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ABSTRACT 


i  .  -y"  • 

^  A  program ,  CORROS,  has  been  written  to  run  on  the  DEC  20  computer  at-DREA.  It  provides  a  means 
'for  analyzing  polarization  resistance  data  from  potentiodynamic  polarization  experiments,  in  order  to  determine 
corrosion  current  densities.  CORROS  accepts  data  files  which  are  transferred  from  an  EG&G  PARC  Model  350 
Corrosion  Measurement  System  to  the  DEC  20  computer  with  another  computer  program,  <Staal>CllA201.  A 
nonlinear  least  squares  curve  fitting  technique  is  used  to  fit  a  curve,  which  satisfies  the  Stem-Geary  equation,  to  the 
experimental  data.  The  corrosion  potential  4^^ ,  the  anodic  and  cathodic  Tafel  constants  ba  and  bc ,  the 

polarization  resistance  Rp,  and  the  corrosion  current  density  i^^  along  with  their  estimated  errors,  and  the 
relative  RMS  error  of  the  fitted  data,  are  determined  and  stored  in  a  file  LIST.  Another  file  PLOTS,  which  is 
accepted  by  <Staal>SAPLTF2  for  graphic  presentation  of  the  data,  is  also  output 
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Le  programme  CORROS  a  et6  6crit  pour  fetre  exploit^  sur  l'ordinateur  DEC  20  du  CRD  A.  Ce  programme 
permet  d'analyser  les  donnees  sur  la  resistance  k  la  polarisation  ob tenues  lors  d'experiences  de  polarisation 
potentiodynamiques,  pour  determiner  ensuite  les  densites  de  courant  de  corrosion.  Le  CORROS  accepte  des 
fichiers  de  donnees  qui  sont  transf6res  d'un  systbme  de  mesure  de  la  corrosion  EG&G  PARC,  module  350,  k 
l'ordinateur  DEC  20  k  l'aide  d'un  autre  infoprogramme,  le  <Staal>Cl  1 A201.  Une  technique  non  lineaire  de  lissage 
par  la  methode  des  moindres  carres  est  udlisee  pour  ajuster  une  courbe,  qui  satisfait  l'equation  de  Stem-Geary,  aux 
donndes  experimentales.  Le  potentiel  de  corrosion  Jes  constantes  de  Tafel  anodique  et  cathodique  ba  et  bc,  la 
resistance  k  la  polarisation  Rp  et  la  densite  de  courant  de  corrosion  i^^,  ainsi  que  les  erreurs  estim6es 
correspondantes  et  l'erreur  quadratique  moyenne  des  donnees  ajustees,  sont  determines  et  mis  en  memoire  dans  le 
fichier  LIST.  Un  autre  fichier,  le  PLOTS,  accepte  par  le  programme  <Staal>SAPLTF2  pour  la  representation 
graphique  des  donnees,  est  aussi  produit 
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Corrosion  rates  can  be  determined  from  weight  loss  measurements3  or,  more  rapidly,  by  less  direct 
electrochemical  techniques.  These  latter  techniques  are  based  on  a  classic  paper  which  correlates  the  use  of  anodic 
and  cathodic  polarization  data  and  mixed  electrode  reactions  to  the  steady  state  (corrosion)  potential  and  the  steady 
state  (corrosion)  current  density4.  The  use  of  electrochemical  techniques  few  determining  corrosion  rates  and 
conducting  corrosion  studies  has  been  rigorously  debated5,  rejected6*7,  and  praised8*9. 

Electrochemical  techniques  are  characterized  by  the  potential  range  over  which  the  data  are  measured  and  by 
the  method  of  data  analysis.  Polarization  resistance  measurements  usually  cover  a  small  potential  range  of  50  mV 
around  the  corrosion  potential1®.  The  linear  polarization  method  was  initially  used  for  the  analysis  of  polarization 
resistance  data  in  order  to  determine  the  slope  of  a  supposed  linear  portion  of  the  polarization  curve11*12.  The 
value  of  this  slope  has  been  called  the  polarization  resistance  (Rp)  and  is  inversely  proportional  to  both  the 

corrosion  current  density  and  the  corrosion  rate.  This  method  has  more  recently  been  strongly  criticized13'16  and 
has  resulted  in  the  development  of  alternative  analysis  techniques. 

One  substitutional  method  which  was  introduced  allowed  simultaneous  calculation  of  the  polarization 
resistance,  the  Tafel  slopes,  the  corrosion  potential  and  the  corrosion  current  density  by  graphic  analysis17  and 
later  by  computer  aided  methods18'19.  A  variety  of  computer  programs  have  been  written  using  linear19  and 
nonlinear20  least  squares  curve  fitting  techniques,  a  simplified  three  point  method21*22,  and  a  four  point 
method23.  Programs  designed  for  use  in  industrial  corrosion  monitoring  systems  and  which  therefore  do  not 
require  sophisticated  computer  analysis  have  also  been  developed24.  A  program  was  written  for  the  analysis  of 
curves  which  deviate  substantially  from  the  family  of  curves  allowed  by  the  Stem-Geary  equation  because  of 
concentration  polarization25.  Another  was  written  specifically  for  microcomputers26  while  yet  another  was  based 
on  a  visual  superposition  of  empirical  and  calculated  curves  through  an  interactive  computer  graphics 
procedure27*28.  Each  of  these  programs  has  unique  approaches,  advantages  and  pitfalls.  Very  few  unbiased 
comparisons  of  the  different  computer  programs  have  taken  place.  It  has  been  demonstrated,  however,  that  the 
least  squares  method  of  analysis  is  more  accurate  than  the  point  methods  of  calculation29. 

The  computer  program  CORROS  was  developed  at  DREA  to  run  on  the  DEC  20  computer  in  order  to 
specifically  analyze  polarization  resistance  data  obtained  through  the  EG&G  PARC  Model  350  Corrosion 
Measurement  System.  The  EG&G  PARC  Model  350  Corrosion  Measurement  System  has  its  own  software  and  is 
capable  of  determining  the  polarization  resistance  Rp  from  polarization  resistance  data  but  does  so  by  the  linear 
polarization  method.  CORROS,  on  the  other  hand,  is  capable  of  simultaneously  determining  the  corrosion 
potential,  the  Tafel  constants,  the  polarization  resistance  and  the  corrosion  current  density  from  polarization 
resistance  data  obtained  from  the  mixed  potential  region  around  the  corrosion  potential.  The  method  is  based  on  a 
nonlinear  least  squares  curve  fit  of  the  data  to  the  Stem-Geary  equation.  The  output  includes  the  corrosion 
parameters,  their  standard  deviations,  and  the  relative  root  mean  square  error  of  the  fitted  data. 


2.  A  DESCRIPTION  OF  CORROS 

CORROS,  as  listed  in  Appendix  A,  is  written  in  the  Pascal  computer  language  and  can  be  easily  divided 
into  three  sections:  dataheader,  datagenerator  and  datafit 

Dataheader  is  designed  to  read  and  recognize  certain  features  in  a  data  file.  The  data  is  transferred  from  the 
EG&G  PARC  Model  350  Corrosion  Measurement  System  to  the  DREA  DEC  20  computer  through  the  use  of 
<Staal>CllA201,  as  listed  in  Appendix  B.  A  data  file,  as  listed  in  Appendix  C,  is  made  up  of  two  parts.  The 
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first  part  contains  information  about  the  experiment,  such  as  the  sample  number,  the  date,  the  tape  and  file 
numbers  where  the  data  is  stored,  the  surface  area,  the  potential  range  of  the  data,  the  potential  scan  rate,  the 
initially  determined  corrosion  potential,  a  value  for  the  polarization  resistance  calculated  by  the  EG&G  PARC 
Model  350  Corrosion  Measurement  System,  and  the  final  determination  of  the  corrosion  potential.  The  second 
part  of  a  data  file  is  the  actual  experimental  results:  pairs  of  potential/current  density  measurements.  Dataheader 
recognizes  the  end  of  the  first  part  of  a  data  file  and  the  beginning  of  the  second  part. 

The  second  section  of  CORROS,  datagenerator,  then  takes  over.  It  reads  the  data  from  a  file,  converts  the 
potential  measurements  from  volts  to  millivolts  and  the  current  measurements  from  nanoamps  per  square 
centimeter  to  milliamps  per  square  centimeter,  and  then  passes  a  cubic  spline  through  the  data  to  allow  access  to 
points  between  the  real  data  points. 

For  a  typical  polarization  resistance  experiment  with  a  potential  range  of  50  millivolts  (<t>corT  ±  25  mV), 

there  are  99  data  points.  Each  data  point  is  a  potential  value  in  volts  and  a  current  density  value  in  nanoamps  per 
square  centimeter.  The  potential  value  of  each  data  point  increases  by  0.5  millivolts,  although  the  EG&G  PARC 
Model  350  Corrosion  Measurement  System  does  not  report  the  points  in  this  manner.  It  instead  rounds  off  the 
potential  values  to  the  nearest  millivolt  or  0.001  volts.  Therefore,  it  gives,  what  appears  at  first  glance  to  be,  two 
pairs  of  data  points  increasing  in  potential  value  by  1  millivolt. 

Datagenerator  interprets  the  data  file  correctly  so  that  each  data  point  increases  in  potential  value  by  0.5 
millivolts.  It  compares  the  potential  value  of  each  data  point  with  the  potential  value  of  the  data  point  read 
previously.  If  the  two  potentials  are  equal,  it  adds  a  correction  of  0.0005  volts  to  the  potential  value  of  the  data 
point  read  most  recently.  For  example,  it  would  interpret  the  following  data  points. 
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Datagenerator  then  stores  the  data  in  arrays  named  "potential"  and  "current"  with  units  of  millivolts  and 
milliamps  per  square  centimeter.  The  ascension  of  the  data  is  checked  to  ensure  that  there  are  no  mistakes  which 
occasionally  appear  while  the  data  is  being  transferred  from  the  EG&G  PARC  Model  350  Corrosion  Measurement 
System  to  the  DEC  20  computer. 

The  data  is  then  interpolated  with  a  cubic  spline  through  an  external  Fortran  subroutine,  SPLINE30.  This 
subroutine  is  listed  in  Appendix  D.  The  purpose  of  interpolation  in  this  program  is  to  have  a  fast  algorithm  to 
calculate  current  density  values  for  potential  values  not  in  the  table  of  data.  The  program  requires  interpolated 
values  of  potential  and  current  density  between  the  real  data  points  at  a  later  part  of  the  program. 

The  subroutine,  SPLINE,  uses  a  technique  called  piecewise  polynomial  intetpolation31.  For  a  set  of  data 
points  or  knots,  X(N)  and  Y(N),  where  N  is  the  number  of  data  points,  X  are  the  abscissas  of  the  data  points,  and 
Y  are  the  ordinates  of  the  data  points,  the  coefficients  B(I),  C(I),  and  D(I),  where  I=1,2,3,...,N,  are  computed  for  a 
cubic  interpolating  spline 

$(X)  *  YCl)  +  B(I)[X-X(I)]  +  C(I)[X-X(I)]2+  D(I)[X-X(I)]3  (1) 
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for 


X(l)  <  X  <  X(l  +  1) 


(2) 


The  third  section  of  CORROS,  datafit,  first  determines  the  corrosion  potential  and  then  screens  the  data  to  a 
symmetric  range  around  this  potential.  The  data  range  around  the  initial  value  of  the  corrosion  potential, 
determined  by  the  EG&G  PARC  350  Corrosion  Measurement  System,  is  25  millivolts  on  either  side.  The 
corrosion  potential  value  actually  shifts  though,  during  the  polarization  of  the  sample.  It  is  consistently  as  much 
as  7  millivolts  cathodic  to  the  initially  determined  value  of  the  corrosion  potential.  The  final  data  set  is  therefore 
lopsided  with  respect  to  the  corrosion  potential.  The  cathodic  range  is  typically  18  millivolts,  while  the  anodic 
range  is  typically  32  millivolts. 


The  algorithm  in  the  program  CORROS  is  written  to  handle  only  symmetric  sets  of  data;  i.e.,  equal  ranges 
of  data  on  either  side  of  the  corrosion  potential.  Only  a  portion  of  the  data  set  is  therefore  analyzed,  typically  18 
millivolts  on  either  side  of  the  corrosion  potential.  This  means  that  up  to  30%  of  the  data  from  an  experiment  is 
ignored. 

The  program  uses  two  external  Fortran  subroutines,  EVAL  and  DEVAL  (from  SEVAL30),  also  listed  in 
Appendix  D,  to  find  where  the  rate  of  change  of  potential  with  current  on  the  spline  is  a  minimum.  The  current 
changes,  at  this  point,  from  travelling  in  a  cathodic  direction  to  an  anodic  direction.  EVAL  is  used  to  evaluate  the 
spline  and  DEVAL  is  used  to  determine  the  spline's  derivative  at  a  given  ordinate  value. 

The  program  continues  with  datafit  by  determining  an  initial  approximation  for  the  Tafel  constants  ba  and 
bc,  the  corrosion  current  density  icorp  and  the  polarization  resistance  Rp,  and  then  by  improving  this 
approximation  to  arrive  at  the  best  solution. 


CORROS  determines  its  initial  approximations  of  the  desired  parameters  through  the  simplification  of  the 
Stem-Geary  equation 
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j  _  exp  ( u)  -  exp (-A»  v) 
Rp(u  +  v) 


.  exp  ( A$  a)  sinh(  A$  v) 
l=  Rpv 


COfr  In  1 0  R„ 


(h*0“  rd(u+ 


v)  2  RpW 


Then,  combining  equations  (10)  and  (11)  gives 


•  =  2  i  exp  (A$a)  3inh(A$v) 


The  properties  of  symmetry  and  asymmetry  are  then  used  on  equation  (12)  to  give 
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which  contains  only  the  unknown  value  of  a. 


The  determination  of  a  at  various  values  of  A<{>  should  give  an  indication  of  the  value  of  A<|>  at  which  the 
Tafel  constants  become  stable.  As  indicated  above,  the  Tafel  constants  are  a  function  of  the  value  of  a,  as  a  result 
of  substitutions  (5),  (6),  (7),  and  (8). 

The  program  is  written  to  search  for  the  minimum  value  of  A<p  at  which  a  becomes  stable,  starting  at  5 
millivolts.  The  value  of  a  does  not  normally  stabilize  until  A<)>  is  typically  8  millivolts.  If  a  does  not  approach  a 
constant,  the  data  are  suspect. 

The  properties  of  symmetry  and  asymmetry  are  again  used  on  equation  (12)  to  give 


[i(A»).-i(-A»)f2.  3lnh'A*v) 


which  contains  only  the  parameters  RD  and  w.  Also, 


(15) 


[i  ( A ♦  )  ■  -i(-A$)]  =  2  icorr  sinh(A^w) 

which  contains  only  the  parameters  and  w. 

In  order  to  separate  out  the  w  value,  the  previously  calculated  spline  is  used  to  evaluate  l/2A<t>  and  3/2A<t>. 
This,  along  with  the  identity 
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where 


sinh  ( A$  w) 


Equation  (17)  is  a  quadratic  in  the  variable 
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and  the  most  positive  root  is  taken  to  calculate  values  of  w. 

The  program  calculates  values  of  w  for  the  data  above  the  minimum  value  of  A<t>  at  which  the  a  values 
become  constant  The  scatter  of  the  w  values  also  gives  an  indication  of  the  quality  of  the  data.  The  average  value 
of  the  w  values  found  is  used  in  equations  (14)  and  (15). 

Values  of  Rp  are  generated  as  a  function  of  Ad>  using  equation  (14).  The  average  values  of  a,  w  and  Rp  are 
then  used  to  generate  initial  approximations  of  ba,  bc,  and  i(.0IT  where 


(22) 


corr  2  Rp  w 

These  approximations  are  then  used  as  the  initial  parameters  for  a  nonlinear  least  squares  fitting  routine 
which  calculates  the  gradient  and  the  Jacobian  matrix  to  the  sum  of  squared  errors,  generates  a  correction  to  each 
parameter,  and  returns  improved  values.  Nonlinear  optimization3^  by  least  squares  is  based  on  the  formulation  of 
the  sum  of  the  squared  errors  as  some  function.  The  gradient  of  this  function  (simply  the  derivative  in  one 
dimension)  is  then  set  equal  to  zero. 

When  the  error  function  is  a  linear  combination  of  the  constants  sought  (4^^  i^  ba,  bc  and  Rp), 

setting  the  gradient  equal  to  zero  results  in  a  set  of  linear  equations  to  be  solved.  When  the  error  function  is  a 
nonlinear  equation  (such  as  for  the  algorithm  in  CORROS),  the  set  of  equations  to  be  solved  are  also  nonlinear. 
An  iterative  solution  is  therefore  needed. 

Newton's  method33  is  one  means  of  arriving  at  this  solution.  In  this  method 


*n+l  "  rn  - 


JAL 

f'('n) 


where  rj  represent  approximations  to  the  root  of  f(x)=0.  Each  new  iterate  rn+ j  is  found  as  the  unique  zero  of  the 
tangent  line  to  the  curve  y=f(x)  at  rn. 

Newton's  method  can  be  expanded  into  more  than  one  dimension.  The  problem,  then,  is  to  find  the  zero  of 
the  gradient.  In  more  than  one  dimension,  the  gradient  is  a  vector  and  its  derivative  becomes  a  matrix,  called  the 
Jacobian.  The  solution  is  also  a  vector. 

The  iteration  can  be  symbolized  as 

Cl  =  $  '  ^  \  (24) 

where  r,  represents  approximations  to  the  solution,  gj  represents  the  gradient  at  i,  and  Jj'1  represents  the  inverse 
Jacobian  evaluated  at  i. 


Simplification  gives 


J  •  A  f  «  -  b; 
n  n 


.  -j 

Ar  =  r. 


rn+1  '  'n 


In  other  words,  Ar  is  the  correction  that  must  be  applied  to  rn  to  generate  rn+] . 


The  problem  is  reduced  to  generating  the  gradient  and  the  Jacobian  and  choosing  the  variables  with  which 
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to  work.  The  following  familiar  equation  is  used. 


In  1 0  A  $ 


■  In  1 0  A  4> 


A  $  =  $  -  $c 


This  equation  contains  4>corr,  ba,  bc  and  i^,^  Rp  is  found  from 


lnl0lcw[  ’c  l]"  K 


Since  there  are  orders  of  magnitude  difference  between  the  actual  parameters,  ba,  bc,  and  i^^,  this 
routine  accepts  these  values  as  constants  and  creates  its  own  set  of  variables,  epsilon  (£),  alpha  (OC),  beta  (P),  and 
gamma  (y).  These  new  variables  are  put  into  the  equation  for  i  as  multipliers  of  <>coir  ba,  bc,  and  icorr 
respectively.  This  results  in 


*  =  V  iCorr 


In  1 0  ( $  - 1 

;*p — =5 


•U  ....  -mio (♦  -c»COT) 

;  exp 


They  are  choosen  this  way  so  that  each  of  their  initial  values  will  be  1  and  any  size  effect  that  would  have  existed 
if  <t>corr,  ba,  bc,  and  ic0rr  were  used  directly  is  negated.  After  each  successful  iteration,  the  constants  <)>coir  ba,  bc, 
and  imn  are  corrected  by  multiplication  by  their  corresponding  new  variable  and  the  new  variables  each  assume  a 

value  of  1  for  the  start  of  the  next  iteration.  This  procedure  is  continued  until  E,  a,  P,  and  y  remain  constant  at  a 
value  of  1  or  the  routine  reaches  its  maximum  allowed  number  of  iterations. 

The  gradient  and  Jacobian  are  calculated,  during  each  iteration,  in  terms  of  £,  a,  P,  and  y.  Since  these 
correction  multipliers  remain  essentially  at  a  value  of  1  throughout,  their  values  do  not  figure  greatly  in  the 
formulae. 

The  gradient  is  calculated  as  the  sum  of  the  partial  derivatives  for  each  correction  multiplier  over  the  data  in 
the  range  being  studied.  The  Jacobian  is  approximated  by  divided  differences  and  is  also  symmetrized  by  averaging 
off  diagonal  terms. 

Two  external  Fortran  subroutines,  DECOMP  and  SOLVE34,  which  are  listed  in  Appendix  D,  are  used  to 
solve  the  matrix  equation  for  the  correction  vector.  Together  they  perform  an  LU  decomposition  on  a  matrix  and 
return  the  condition  number  of  the  matrix.  The  coefficient  matrix  is  passed  to  DECOMP,  and  SOLVE  uses  the 
LU  form  along  with  the  solution  vector  to  find  a  solution.  The  correction  vector  is  restricted  by  two  conditions: 


i)  £,  a,  (3,  and  y  are  not  allowed  to  change  by  more  than  50%  per  iteration,  and 

ii)  the  root  mean  square  error  at  the  corrected  estimate  must  be  less  than  the  root  mean  square 
error  at  the  start 

This  protects  the  routine  from  running  away  from  or  overshooting  an  error. 

The  desired  constants  are  then  output  with  their  estimated  error  and  the  predicted  data  are  tabulated  with  the 
actual  data.  The  software  can  be  easily  modified  so  that  the  data  is  passed  to  another  file  which  is  suitable  for 
graphic  presentation,  or  to  a  file  which  just  contains  results. 


3.  RUNNING  CORROS 

There  are  two  versions  of  a  CORROS  run.  An  example  of  the  simplified  version  is  in  Appendix  E  and  the 
detailed  version  is  in  Appendix  F.  The  difference  between  the  two  versions  is  the  amount  of  output  to  the 
terminal.  During  a  detailed  run,  the  program  outputs  the  gradient,  the  Jacobian  and  its  construction,  the  correction 
vector,  and  its  treatment  at  each  iteration.  During  a  simplified  run,  the  program  operates  the  same  routines  but 
suppresses  the  output  The  results  and  the  information  supplied  to  PLOTS  and  LIST  are  identical. 

The  program  is  initiated  by  ex@corrode  which  results  in  the  loading  of  the  program  CORROS  and  the 
subroutines  SPLINE30,  SEVAL30,  DECOMP34  and  SOLVE34,  which  are  found  in  CORSUB.FOR. 

@ex@corrode 

PASCAL:  CORROS 

FORTRAN:  CORSUB 

SPUNE 

SEVAL 

DECOMP 

SOLVE 

LINK:  Loading 

l LNKXCT  CORROS  execution ] 

Execution  starts  with 


INPUT  : 

after  which,  the  user  responds  with  a  carriage  return,  then 

OUTPUT  : 

after  which,  the  user  again  responds  with  a  carriage  return,  and  then 

PLOTS  : 


The  user  now  has  a  choice  of  either  naming  the  output  file,  suitable  for  graphic  presentation,  PLOTS  by 
responding  with  a  carriage  return  or  giving  this  file  another  name  such  as  P05F04.DAT  for  example.  A  listing  of 
a  PLOTS  type  file  is  found  in  Appendix  G.  The  next  output  is 


DATA  : 


The  user  must  now  give  the  name  of  the  data  file,  such  as  T05F04DAT,  that  is  to  be  analyzed.  After  the  next 
output 

UST  : 

the  user  again  has  a  choice  of  either  naming  the  output  file,  which  contains  a  summary  of  results,  LIST  by 
responding  with  a  carriage  return  or  giving  this  file  another  name  such  as  L05.DAT  for  example.  A  listing  of  a 
LIST  type  file  is  found  in  Appendix  H. 

The  program  next  asks  the  user  if  they  would  like  to  read  the  introduction. 

Would  you  like  to  read  the  introduction  ?  yes  =  1 

If  the  user  wishes  to  do  so,  they  respond  with  a  1. 

1 

A  program  to  find  the  corrosion  potential  Ecorr,  the  anodic  and  cathodic  Terfel  constants  ba  and  be,  the  polarization 
resistance  Rp,  and  the  corrosion  current  density  Icorr.from  the  mixed  potential  region  around  Ecorr.  The  data  is 
read  from  polarization  resistance  data  from  the  EG&G  PARC  Model  350  Corrosion  Measurement  System  by 
another  program,  <Staal>Cl  lA20f  The  program  can  be  run  in  two  modes:  detailed  and  simple.  The  only 
difference  being  the  amount  of  output  to  the  terminal.  In  detailed  mode  the  program  outputs  the  gradient,  the 
Jacobian  and  its  construction,  the  correction  vector  and  its  treatment  at  each  iteration.  The  simple  mode  operates 
the  same  routines  but  suppresses  output. 

If  not,  the  response  is  anything  except  1.  The  program  then  asks  the  user  if  a  detailed  run  is  required. 

do  you  want  a  detailed  run  ?  yes  =  1 
Again  1  means  yes  and  anything  else  except  1  means  no. 

After  the  program  reports  the  information  which  is  in  the  data  header  at  the  beginning  of  the  data  file, 


SAMPLE  4  DATE  28.10  TAPE  2  FILE 

MV/SEC  0.167  ECORR  -0.252 

RESULTS  RP  5.257E3 

V’.V, 

ECORR -0.256  ***  end  of  tape  header  *** 

and  the  number  of  data  points  in  the  data  file, 

• » ■ . 

99  data  points  read  in 

it  asks  the  user  if  they  would  like  to  see  a  tabulation  of  the  data. 

do  you  want  to  see  the  data  tabulated?  yes  =  1 
1  again  means  yes  and  anything  else  means  no. 


m 

..\v; v; 

r*  VVV a  kV  A  .  »  - 


>•  ■  „•  .*  /  -■  /  V  V  V  ••  .4**  ‘  „•  * 


*  .  ‘  •  *  •  *  •  ’  *  *  *  •  *  a  *  i*  -  '  • 


The  program  reports  the  initial  guess  for  4^^, 
estimated  E  '”r  is  -255.11 
the  size  of  one  half  of  the  symmetric  range  around 

range  around  Ecorr  reduced  to  209 

and  the  number  of  data  points  in  the  symmetric  range  which  are  being  used  in  the  analysis. 

84  points  in  range 

The  program  next  reports  the  value  of  a  at  each  0.5  mV  interval  of  A4>,  starting  at  A4>  =  5.0  and  going  to  the  end  of 
the  symmetric  range  of  data. 


atoP  =  5.0 
at  oP  =  6.0 

a=  -0.000173 
a  =  0.002409 

cuoP  m  5.5 
at  oP  =  65 

a  «  0.001862 
a  =  0.003409 

at  oP  =  19.0 
at  oP  =  20.0 

a  =  0.003753 
a  m  0.001458 

atoP  =  195 
at  oP  =  205 

a  =  0.002325 
a  =  -0.000052 

After  reporting  the  average  value  of  a,  the  program  asks  the  user  for  the  value  of  A4>  at  which  the  user  wants  the 
program  to  continue  with  the  analysis. 

average  a  is  0.006881 

enter  starting  oP:  minimum  oP  where  a  becomes  reasonably  constant 

The  user  normally  chooses  A<>  =  5  but  by  scanning  the  list  of  a  values,  may  also  choose  the  minimum  value  of  A<p 
at  which  the  value  of  a  becomes  reasonably  constant.  This  choice  can  sometimes  solve  nonconvergence  problems 
with  certain  data  files. 

The  program  then  reports  the  average  value  of  w, 
average  w  is  0.057436 
an  improved  reestimate  of  the  value  of  a, 

re-estimate  of  a  is  0.0069 


an  average  value  ■  f  Rp, 

average  Rp  is  521 7. 


an  estimate  of 

I 
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estimated  Icon  is  0.00166878 


l:  *<  *  • 
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a  reestimate  of  the  value  of  w, 

re-estimate  of  w  is  0.05757037 
the  initial  estimates  of  4>corr  ba,  bc,  and  Rp, 

initial  EcorrJt>aJbcJcorr  and  Rp  are 
-255.11  35.73  45.43  0.001665  5217. 

and  the  initial  relative  root  mean  square  error  of  the  fit. 

initial  relative  RMS  error  is  0.082915 

If  the  initial  values  are  acceptable,  the  user  responds  with  a  1. 

enter  1  if  acceptable 

If  they  are  not  acceptable,  the  user  responds  with  anything  else  except  a  1  and  is  then  asked  to  supply  initial 
estimates  for  ba,  bc  and  Rp. 

2 

enter  babe  Rp*  in  mvolts/decade  and  ohms  * 

36455200 

relative  RMS  error  =  0.079201 
enter  1  if  acceptable 
1 

The  program  then  asks  the  user  to  respond  with  a  1  if  they  wish  the  program  to  continue  and  the  iterations  to 
begin,  or  to  respond  with  a  99  if  they  wish  the  program  to  stop. 

enter  1  to  continue  :  99  to  quit 

After  each  iteration,  if  the  run  is  detailed,  the  program  reports  the  gradient,  the  Jacobian  matrix  and  the 
correction  vector. 


gradient  at  each  parameter  +  del 

5559737E-04  -2.110495E-05  2.477400E-05  -6.371798E-06 
5.498207E-04  -2.090320E-05  2.471899E-05  -6.459422E-06 
5500678E-04  -2.091288E-05  2.472174E-05  -6.454742E-06 
5500761E-04  -2.092529E-05  2.471082E-05  -6.436885E-06 

gradient:  dEcorr  dbadbcdl con- 

55  00043  E-04  -2.091514E-05  2.471674E-05  -6.449062E-06 
Jacobian  before  symmetry  check 

5.969370E-01  -1.83681 5E-02  6.341 725 E-03  7.173367E-03 
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-1.898122E-02 

5.726088E-03 

7.726374E-03 

Jacobian  matrix 
5.969370E-01 
-1.867469E-02 
6.033906E-03 
7.449870E-03 


1.193507E-03 

2.246452E-04 

-1.036022E-03 


-1.867469E-02 
1.193507E-03 
2. 249294 E-04 
-1.02S759E-03 


2.252136E-04 

4.996309E-04 

-5.680192E-04 


6.033906E-03 

2.249294E-04 

4.996309E-04 

-S.801297E-04 


-1.015496E-03 

-5.922402E-04 

1.217643E-03 


7.449870E-03 

-1.025759E-03 

-5.801297E-04 

1.217643E-03 


3246053E-01 


initial  correction  vector 

5.176613E-03  3579467E-01  1.037736E-01 
error  overshoot:  reduced  correction  by  half 


The  program  then  reports  improved  estimates  of  4>corr.  ba,  bc,  and  Rp,  and  the  reduced  relative  root  mean 

square  error  of  the  fit  for  both  detailed  and  simplified  runs.  The  user  may  stop  the  iterations  at  any  time  by 
responding  with  a  99  or  may  continue  by  responding  with  a  1 . 


Ecorr  ba  bc  Icorr  Rp 
-255.77  42.12  47.78  0.001935  5024. 

relative  RMS  error  now  =  0.073427 
enter  1  to  continue  :  99  to  quit 
1 


The  program  reports  when  there  is  no  further  improvement  available,  the  user  responds  with  a  99,  and  the  program 
reports  the  final  relative  root  mean  square  error  of  the  fit,  the  final  estimates  of  ba,  bc,  icorT  and  Rp,  and  the 
error  in  each  estimated  parameter. 

no  further  improvement  available 

Ecorr  ba  bc  Icorr  Rp 
-256.18  34.67  33.74  0.001368  5430. 

relative  RMS  error  now  =  0.060409 
enter  1  to  continue  :  99  to  quit 
99 

relative  RMS  error  now  =  0.060409 

Ecorr  «  -256.18  millivolts 

ba  =  34.67  +!-  4.640  millivolts 

bc  =  33.74  +1-  2.949  millivolts 

Rp  =  5430.  +/-  460.8  ohms 

Icorr  =  0.001368  +/-  0.0001 1417  millamps  I  err? 

The  program  finally  asks  the  user  for  the  tape  number  of  the  data  file  and  the  file  number  of  the  data  file. 

Referring  to  the  example,  the  user  would  respond  with  5  and  4  respectively. 

TAPE  number  ? 

5 

FILE  number  ? 

4 


A  summary  of  the  results  of  the  analysis  is  found  in  LIST  or  under  a  file  name  that  the  user  has  chosen. 
This  file  is  appended  so  it  may  contain  summaries  of  previous  analysis  results.  The  data  and  the  fitted  curve  to  the 
data  may  be  graphically  presented  together  by  running  <Staal>SAPLTF2  and  by  giving  this  program  the  name  of 
the  file  which  contains  the  analysis  suitable  for  graphic  presentation;  either  PLOTS  or  a  name  chosen  by  the  user. 
The  results  of  the  analysis  of  data  file  T05F04.DAT  are  shown  in  Figure  1. 


4.  CORROS  VERSUS  TAFEL  PLOT  ANALYSIS 

CORROS  is  also  capable  of  analyzing  polarization  data  over  a  range  larger  than  that  of  polarization 
resistance  data.  Donahue^  presented  an  illustration  of  the  analysis  of  electrochemical  rate  data  for  activation 
controlled  reactions,  using  a  Tafel  plot  technique.  For  the  purpose  of  a  comparison,  the  illustration  was  also 
analyzed  with  CORROS.  Donahue's  data  is  listed  in  Appendix  I.  The  analysis  of  the  data,  using  a  simplified  run 
of  CORROS,  is  listed  in  Appendix  J.  The  output  PLOTS  file  and  the  output  LIST  file  are  listed  in  Appendices  K 
and  L  respectively.  Donahue’s  Tafel  plot  analysis  is  illustrated  in  Figure  2  while  the  CORROS  analysis,  showing 
the  data  and  the  fitted  curve,  is  presented  in  Figure  3. 

The  results  of  the  two  analysis  techniques  are  significantly  close  as  is  shown  in  Table  I. 


5.  CONCLUDING  REMARKS 

The  computer  program  CORROS  provides  a  means  for  the  analysis  of  polarization  resistance  data  from 
electrochemical  corrosion  experiments,  for  the  purpose  of  simultaneously  determining  corrosion  current  densities, 
Tafel  slopes,  polarization  resistance  values  and  corrosion  potentials.  The  results  of  the  analysis  are  summarised  in 
one  output  file  and  are  readily  available  for  graphic  presentation  in  another  output  file.  CORROS  may  be  easily 
used  by  persons  with  no  detailed  knowledge  of  corrosion  science  and/or  no  experience  with  computers. 


TABLE  I 


L  PLOT  ANALYSIS  AND  THE  CORROS  ANALYSIS 
)NAHIJE'S  ILLUSTRATION 


ba(mV)  bc  (mV)  Rp  (k!2)  (mA/cm2) 


40.  120. 


1.0 


40.57±0.71  126.70+1 3. 15  0.01±0.00 


1.07±0.03 


-2/0  -260  -250 

potential  (millivolts  vs. 


-240 

SCE) 


Figure  1  <Staal>SAPLTF^  Presentation  of  the  CORROS  Analysis  of  T05F04.DAT 


Figure  3  CORROS  Analysis  of  Donahue's35  Illustration 


APPENDIX  A 


CQRRQS 


PROGRAM  CorrosionAnalysis(INPUT:/, OUTPUT  JJLOTSJ)ATA,LIST); 

{  A  program  to  find  the  corrosion  potential  Ecorr,  the  anodic  and  cathodic  tafel  constants  ba  and  be,  the  } 
{  polarization  resistance  Rp  and  the  corrosion  current  density  Icorr,  from  the  mixed  potential  region  } 
{  around  Ecorr.  The  data  is  read  from  polarization  resistance  data  from  an  EG&G  PARC  Model  350  } 

{  Corrosion  Measurement  System  by  another  program,  <Staal>Cl  1A201.  The  program's  basic  routine  } 
{  is  a  nonlinear  least  squares  fit  of  the  data.  Output  includes  the  corrosion  parameters,  their  estimated  } 
{  errors,  and  the  relative  RMS  error  of  the  fitted  data.  These  parameters  can  be  found  in  a  file,  LIST..  } 
{  The  real  data  and  the  fitted  data  are  also  output  to  another  file,  PLOTS..,  and  a  plot  can  be  output  by  } 
{  running  <STAAL>SAPLTF2.  } 

CONST 

numdata  =  200; 
order  =  4; 


TYPE 

index  =  L.numdata; 

sindex  =  1.. order; 

svector  =  ARRAY[sindex]  of  real; 

isvector  =  ARRAY[sindex]  of  integer; 

vector  =  ARRAY  [index]  of  real; 

matrix  =  ARRAY  [sindex, sindex]  of  real; 


VAR 

DATA.PLOTS.LIST  :  text; 

gradstorejac  :  matrix; 

gradient,mgradient,work  ;  svector; 

potential, current, voltdata,currdata,BB,CC,DD  :  vector; 

ipvt  :  isvector; 

i,j,m,n,numofpts,count,ans,asn  :  integer; 

Rp,Ecorr, Icorr, ba,bc,cond  :  real; 

detailed  ;  boolean; 


PROCEDURE  dataheader; 

CONST 
size  =  6; 
blank  = ' '; 


TYPE 


wordtype  =  ARRAY[l..size]  of  char; 
testtype  =  PACKED  ARRAY[l..size]  of  char; 


ch  :  char; 

i  :  integer; 

testwordjceyword  :  testtype; 
seive  :  wordtype; 

PROCEDURE  skipblanks; 

BEGIN 

REPEAT  read(DATA,ch)  UNTIL  (ch  <>  blank); 
END; 

PROCEDURE  initseive; 

BEGIN 

skipblanks; 
seive[l]  :=  ch; 

FOR  i  :=  2  TO  size  DO 
BEGIN 

read(DATA,ch); 
seive[i]  :=  ch; 

END; 

END; 

PROCEDURE  Iookfor(keyword:testtype); 

BEGIN 

initseive; 

REPEAT 

write(seive[l]); 

read(DATA,ch); 

FOR  i  :=  1  TO  size-1  DO 
seive[i]  :=  seive[i+l]; 
seive(size)  ;=  ch; 
pack(sei  ve,  1  .testword) ; 

UNTIL  (testword  =  keyword); 
writeln; 

FOR  i  :=  1  TO  size  DO 
write(seive[i]); 

END; 

BEGIN  { dataheader } 
writeln; 

('  Would  you  like  to  read  the  introduction  ?  yes  =  1'); 
readln;  read(ans); 

IF  ans  =  1 
THEN 


writelnC  A  program  to  find  the  corrosion  potential  Ecorr,  the  anodic  and  cathodic  tafel  constants  ba  and  be,  ’) 
writelnO  the  polarization  resistance  Rp,  and  the  corrosion  current  density  Icorr,  from  the  mixed  potential  ') 

writelnC  region  around  Ecorr.  The  data  is  read  from  polarization  resistance  data  from  an  EG&G  PARC  ') 

writelnC  Model  350  Corrosion  Measurement  System  by  another  program,  <STAAL>C1 1 A20.  The  ’) 

writelnC  program  can  be  run  in  two  modes:  detailed  and  simple.  The  only  difference  being  the  amount  of  '): 

writelnC  output  to  the  terminal.  In  detailed  mode  the  program  outputs  the  gradient,  the  Jacobian  and  its  ') 
writelnC  construction,  the  correction  vector  and  its  treatment  at  each  iteration.  The  simple  mode  operates  '); 
writelnC  the  same  routines  but  suppresses  output.  '): 


END; 

writeln; 

writelnC  do  you  want  a  detailed  run  ?  yes  =  1 '); 
readln;  read(ans); 

IF  ans  =  1  THEN  detailed  :=  TRUE 
ELSE  detailed  :=  FALSE; 

reset(DATA); 
lookfor('RESULT') ; 
lookfor(ECORR '); 
read(DATA,ch);  write(ch); 

WHILE  NOT  eoln(DATA)  DO 
BEGIN 

read(DATA,ch);  write(ch); 

END; 

writelnC  ***  end  of  tape  header  ***'); 
writeln; 

END;  { dataheader  } 

(These  routines  are  FORTRAN  subprograms  supplied  at  load  time. 

(All  four  can  be  found  in  CORSUB.FOR. 

PROCEDURE  spline(n  :  integer;  xdata.ydata  :  vector;  VAR  b,c,d  :  vector); 

{  fits  an  interpolating  cubic  spline  to  a  set  of  data  of  strictly  increasing  ordinate  (x)  value  } 

extern  FORTRAN; 

PROCEDURE  seval(n  :  integer;  xval :  real;  xdata,ydata,b,c,d  :  vector;  VAR  yval.dyval :  real); 

{  evaluates  the  spline  (y)  and  its  deravitive  (y')  at  a  given  ordinate  (x)  value  } 
extern  FORTRAN ; 

PROCEDURE  decomp(ndim,n  :  integer;  VAR  A  :  matrix;  VAR  cond  :  real;  VAR  ipvt :  isvector; 

VAR  work  :  svector); 


{  performs  LU  decomposition  of  an  n  x  n  matrix  A  and  returns  the  condition  number  of  the  matrix  } 


extern  FORTRAN; 

PROCEDURE  solve(ndim,n  :  integer;  VAR  A  :  matrix;  VAR  B  :svector;  VAR  ipvt :  isvector); 

{  solves  a  set  of  n  simultaneous  equations  in  n  variables.  The  coeffitient  matrix  A  is  passed  to  } 

{  DECOMP  and  SOLVE  and  uses  the  LU  form  along  with  the  solution  vector  B  to  solve  Ax  =  B  } 

extern  FORTRAN ; 

PROCEDURE  DataGenerator; 

{  takes  raw  data  from  a  fileDAT  and  converts  potenials  from  volts  to  mVolts  and  current  densities  from  } 
{  namps  to  mAmps  per  sqcm  and  passes  a  cubic  spline  through  the  data  to  allow  access  to  points  ’between'  } 
{  data  points.  } 

VAR 

voltage, amperage, vstp :  real; 

PROCEDURE  scale;  {  scales  volts  to  mVolts  and  namps  to  mamps  } 

BEGIN 

FOR  i  :=  1  TO  numofpts  DO 
BEGIN 

potential^]  :=  potential^]  *  1E3; 
current[i]  :=  current[i]  /  1E6; 

END; 

END; 

PROCEDURE  reorder;  {  reorders  decreasing  voltage  data  } 

VAR  tempvolt,tempamp  :  real; 
i,nexti  :  integer; 


BEGIN 

writeln('reordering'); 

FOR  i  :=  1  TO  (numofpts  DIV  2)  DO 
BEGIN 

nexti  :=  numofpts  +  1  -  i; 


tempvolt  :=  potential^]; 
potential^]  :=  potential[nexti]; 
potential[nexti]  :=  tempvolt; 


tempamp  :=  current[i]; 
current[i)  :=  currentfnexti]; 
current[nexti]  :=  tempamp; 


END; 


END;  {  reorder  } 


PROCEDURE  viewdata; 


writeln;  writeln(’do  you  want  to  see  the  data  tabulated?  yes  =  1'); 
readln;read(ans); 

IF  ans  =  1  THEN 
BEGIN 

writeln('****************  ********************  ******************♦****')• 
writeln; 

writeln('  Potential  Current  Potential  Current’); 
writeln; 

FOR  j  :=  1  TO  numofpts  DIV  2  DO 
writeln(j:2, potential  [j]:8:2,current(j];15:6, 

(j  +  numofpts  DIV  2):8,potential[j+numofpts  DIV  2]:8:2, 
current[j+numofpts  DIV  2]:16:6 ); 

writelnC***********************************************************')- 

END; 

END;  {  viewdata  } 

BEGIN  {  datagenerator  } 

i>  l; 

read(DATA,voltage,amperage); 

EF  amperage  <  0  THEN  vstp  :=  0.0005  ELSE  vstp  ;=  -0.0005; 

WHILE  NOT  eofp  ATA)  DO 
BEGIN 

potential!!]  ;=  voltage; 
current[i]  :=  amperage; 
read(D  AT  A,  voltage, amperage); 

{ this  section  detects  data  which  have  the  same  voltage  because  the  EG&G  PARC  Model  350  Corrosion  } 
{  Measurement  System  can  only  pass  3  digits  and  either  adds  or  subtracts  half  of  the  last  digit  to  those  } 

{  voltages  depending  on  whether  the  voltages  are  increasing  or  not  } 

IF  voltage  =  potential^]  THEN  voltage  :=  voltage  +  vstp; 

i  :=  i  +  1; 

END;  {  WHILE  } 

numofpts  :=  i  -  1; 
scale; 

writeln(numofpts:3,'  data  points  read  in '); 

IF  potential^]  >  potential[4]  THEN  reorder; 
viewdata; 

spline(numofpts,  potential,  current,  BB,CC,DD); 


END;  {  datagenerator } 


PROCEDURE  screendata( window  :  real); 


{  screens  out  data  that  is  more  than  a  'window'  away  from  Ecorr  } 

BEGIN 

FOR  i  :=  1  TO  numdata  DO 
BEGIN 

voltdata[i]  :=  0;  currdata[i]  :=  0; 

END; 

j>0; 

FOR  i  :=  1  TO  numofpts  DO 
BEGIN 

IF  abs(potentiaI[i]  -  Ecorr)  <=  window 
THEN 

BEGIN 

j  >j  +  i; 

voltdata[j]  :=  potential^]; 

currdata[j]  :=  currentfi]; 

END; 

END; 
count  :=  j; 

writeln(count4,'  points  in  range  '); 

END;  {  screendata  } 

PROCEDURE  datafit; 

{  let  L  =  In(10),  Ecorr  be  the  corrosion  potential  in  mvolts  and  E  -  Ecorr  be  the  overpotential  oP  in  } 

{  millivolts.  Let  I  be  the  total  current  density  in  mamps  and  Icorr  be  the  corrosion  current  density  in  } 

{  mamps.  Then  the  equation  that  governs  I  wrt  P  is:  1 .  I  =  Icorr  *  (exp(L  *  oP/ba)  -  exp(-L  *  oP/bc))  } 

{  where  ba  and  be  are  the  anodic  and  cathodic  Tafel  constants.  Further,  the  slope  of  the  I  vs  oP  curve  at  } 

{  oP  =  0  is  a  constant  @  oP=0,  (dl/doP)  -  1/Rp,  where  Rp  is  the  polarization  resistance  and  @  } 

{  oP=0,  (dl/doP)  =  Icorr  *  ln(10)  *  (1/ba  +  1/bc).  Thus  2.  1/Icorr  =  L  *  Rp  *  (1/ba  +  1/bc).  Now  } 

{  substitute  u  =  L/ba  and  v  -  L/bc  and  2.  into  1.  and  get  } 

{  3.  I  =  (exp(oP  *  u)  -  exp(-oP  *  v))/(Rp  *  (u  +  v)).  Now  substitute  u  =  w  +  a  and  v  =  w  -  a  and  get  } 

{  4.  I  =  exp(P  *  a)  *  sinh(P  *  w)/(Rp  *  w).  The  variables  in  equation  4  can  be  partly  or  wholly  } 

{  separated.  Consider  I  =  I(oP) :  i)  I(oP)/I(-oP)  =  -exp(2  *  oP  *  a)  contains  only  a,  } 

{  ii)  sqrt(I(oP)*-I(-oP))  =  sinh(oP  *  w)/(Rp  *  w)  only  Rp  and  w,  or  substituting  Icorr  into  (ii),  } 

{  iii)  sqrt(I(oP)M(-oP))  =  2  *  Icorr  *  sinh(oP  *  w).  One  variable  which  has  not  been  dealt  with  yet,  Ecorr,  } 
{  is  hidden  in  P  or  oP  and  since  these  appear  as  arguments  to  the  exponential,  it  is  very  necessary  to  } 

{  accurately  determine  Ecorr.  The  initial  approximation  to  Ecorr  is  the  potential  E  at  which  the  current  } 
{  density  1(E)  =  0.  This  E  is  found  by  solving  for  the  root  of  spln(E)  =  0,  where  spin  is  the  cubic  spline  } 

{  which  interpolates  the  data  and  this  is  done  using  Newton's  method,  xl  =  xO  -  y/y'.  } 

VAR 

V,  oP,  oPincr,  minoP,  wO,  suml,  sum2,  sum3,  sum4,  yval,  dyval,  tempi,  temp2,  L,  kp,  delta, 
w,  sumw,  sdevw,  estEcorr,  range,  a,  suma,  sdeva,  den,  z,  IoP,  Iplus,  Iminus,  lp3half,  lm3half,  oPrange, 
rms,  newrms,  epsilon,  alpha,  beta,  gamma,  sdevlcorr,  sdevba,  sdevbc,  sdevRp  :  real; 


i,j,num,rej :  integer; 

FUNCTION  eval(  V :  real) :  real; 

{  returns,  from  the  data  spline,  a  current  value  for  a  given  input  potential  in  mvolts  } 
BEGIN 

seval(numofpts,V, potential, current,BB,CCJ>D,yval,dyval); 
eval  :=  yval; 

END; 

FUNCTION  deval(  V  ;  real) :  real; 

{  returns  the  rate  of  change  of  current  at  a  potential } 

BEGIN 

seval(numofpts,V,potential,current,BB,CC,DD,yval,dyval); 
deval  :=  dyval; 

END; 

FUNCTION  fitcurrent(V, Ecorr, ba,bc, Icorr :  real) :  real; 

BEGIN 

oP  :=  V  -  Ecorr; 

IF  oP  <>  0 
THEN 

fiteurrent  :=  Icorr  *  (exp(L*oP/ba)  -  exp(-L*oP/bc)) 

ELSE 

fiteurrent  :=  0.0; 

END; 

FUNCTION  rmserr(Ecorr,ba,bc,Icorr :  real) :  real; 

VAR  suml,sum2  :  real; 

i  j  :  integer; 

BEGIN 

suml  :=  0;  sum2  :=  0;  j  :=  0; 

FOR  i  :=  1  TO  count  DO 
BEGIN 

oP  :=  voltdata[i]  -  Ecorr; 

IF  ((abs(oP)  >=  minoP) 

AND  (abs(oP)  <=  oPrange  +  1)) 

THEN 

BEGIN 

suml  :=  suml  +  sqr(currdata[i]  -  fitcurrent(oP  +  Ecorr, Ecorr, ba,bc, Icorr)); 
sum2  :=  sum2  +  sqr(currdata[i]); 


END; 

rmserr  >  sqrt(suml/sum2); 
END; 


PROCEDURE  improve(VAR  Ecorr, ba,bc, Icorr :  real); 


{  This  routine  takes  initial  values  of  Ecorr,  ba,  be,  and  Icoit  and  returns  improved  values,  least  squares  } 

{  sense,  by  applying  a  nonlinear  Newton  iteration  in  four  variables.  The  iteration  is  of  the  form  } 

{ JO  *  dX  =  -G(X)  where  dX  is  the  correction  vector ,  JO  is  the  Jacobian  matrix  evaluated  at  XO  and  GO  } 
{  is  the  gradient  of  the  sum  of  squared  errors  at  XO,  the  initial  parameters.  Since  there  are  orders  of  } 

{  magnitude  difference  between  the  actual  parameters  Ecorr,  ba,  be,  and  Icorr  this  routine  accepts  these  } 

{  values  as  constants  and  creates  its  own  set  of  variables  epsilon,  alpha,  beta,  and  gamma.  These  new  } 

{  variables  are  put  into  the  equation  for  I  as  multipliers  of  Ecorr,  ba,  be,  and  Icorr  respectively.  They  are  } 

{  chosen  this  way  so  that  each  of  their  initial  values  will  be  1  (one)  and  this  negates  any  size  effect  that  } 

{  would  have  existed  if  Ecorr,  ba,  be,  and  Icorr  were  used  directly.  After  each  successful  iteration,  the  } 

{  'constants'  Ecorr,  ba,  be,  and  Icorr  are  corrected  by  multiplication  by  their  corresponding  new  variable  } 

{  and  the  new  variables  each  assume  a  value  of  1  (one)  for  the  start  of  the  next  iteration.  The  procedure  } 

{  is  continued  until  epsilon,  alpha,  beta,  and  gamma  remain  constant  at  a  value  of  1  (one)  or  the  routine  } 

{  iterates  past  its  maximum  allowed  runs  (user  set ,  default  of  50).  } 


VAR 

tempi,  temp2,  predcurr,  error,  pepsilon,  palpha,  pbeta,  pgamma  :  real; 
i,j,m,n  :  integer; 

PROCEDURE  grad(eps,alp,beta,gam :  real); 

VAR  maxgrad :  real; 

BEGIN 

{  'p'  is  used  to  represent  differentiation  } 

pepsilon  :=0;  pgamma  :=  0;  palpha  :=  0;  pbeta  :=  0; 

FOR  i  :=  1  TO  count  DO 
BEGIN 

oP  :=  voltdata[i]  -  eps*Ecorr; 

IF  ((abs(oP)  >=  minoP) 

AND  (abs(oP)  <=  oPrange  +1)) 

THEN 

BEGIN 

predcurr  :=  fitcurrent(voltdata[i],eps*Econ,alp*ba,beta*bc,gam*IcorT); 
error  :=  currdata[i]  -  predcurr; 

tempi  :=  exp(  L  *  oP/ba)/ba; 
temp2  :=  exp(-L  *  oP/bc)/bc; 

pepsilon  :=  pepsilon  +  (2  *  error  *  Icorr  *  L  *  Ecorr  *  (tempi  +  temp2)); 
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pgaxnma  :=  pgamma  +  (-2  *  error  *  predcurr); 

palpha  :=  palpha  +  (2  *  error  *  Icorr  *  L  *  oP  *  tempi); 

pbeta  :=  pbeta  +  (2  *  error  *  Icon  *  L  *  oP  *  temp2); 

END; 

END; 

gradient[l]  :=  pepsilon; 
gradient[2]  :=  palpha; 
gradient[3]  :=  pbeta; 
gradient[4]  :=  pgamma; 

END;  {  grad  } 

PROCEDURE  jacobian(Ecorr,ba,bc,Icorr :  real); 

{  The  Jacobian  matrix  is  made  up  of  the  derivatives  of  the  gradient  in  the  form:  Jac[i,j]  =  dgrad[i]/dX[i]. 
{  This  routine  approximates  the  derivatives  by  a  divided  difference  and  further  assures  that  the  matrix  is 
{  exactly  symmetric,  accounting  for  some  roundoff. 

CONST  del  =  IE-5; 

BEGIN 

grad(l+del, 1,1,1); 

FOR  m  >  1  TO  4  DO 

gradstore[l,m]  :=  gradient[m]; 

grad(  1,1  +del,  1,1); 

FOR  m  :=  1  TO  4  DO 

gradstore[2,m]  >  gradient[m]; 

grad(  1,1,1 -fdel,  1); 

FOR  m  :=  1  TO  4  DO 

gradstore[3,m]  :=  gradient[m]; 

grad(l,l,l,l+del); 

FOR  m  :=  1  TO  4  DO 

gradstore[4,m]  :=  gradient[m]; 

write  In; 

{  Writes  out  the  gradient  at  the  4  points  around  current  estimate,  for  use  in  finding  Jacobian. 

IF  detailed  THEN 
BEGIN 


writeln('  gradient  at  each  parameter  +  del'); 
FOR  m  >  1  TO  4  DO 
BEGIN 

FOR  n  :=  1  TO  4  DO 
write(gradstore[m,n]); 
write  In; 

END; 

writeln; 

END; 

grad(l, 1,1.1); 

IF  detailed  THEN 
BEGIN 

writelnO  gradient:  dEcorr,dba,dbc,dIcorT’); 
FOR  m:=  1  TO  4  DO 
write(gradient(m]); 
writeln;writeln; 

END; 

FOR  m  :=  1  TO  4  DO 
FOR  n  :=  1  TO4DO 

Jat[m,n]  :=  (gradstore[n,m]  -  gradient[m])/del; 

{  Writes  out  the  Jacobian  before  forced  symmetry.  } 

IF  detailed  THEN 
BEGIN 

writelnC  Jacobian  before  symmetry  check '); 
FOR  m  :=  1  TO  4  DO 
BEGIN 

FOR  n  :=  1  TO  4  DO 

write(Jac[m,n]); 

writeln; 

END; 

END; 


writeln; 

FOR  m  :=  1  TO  4  DO 
FOR  n  :=  1  TO  4  DO 
BEGIN 
IF  m  <>  n 
THEN 
IF  m  <  n 

THEN  Jac[m,n]  ;=  (Jac[m,n]  +  Jac[n,m])/2; 

END; 

FOR  m  :=  1  TO  4  DO 
FOR  n  :=  1  TO  4  DO 
IF  m  >  n 

TNI  N  Jac[m,n]  ;=  Jac[n,m]; 
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IF  detailed  THEN 
BEGIN 

writelnC  Jacobian  matrix '); 

FOR  m:=  1  TO  4  DO 
BEGIN 

FOR  n  :=  1  TO  4  DO 

write(Jac[m,n]); 

writeln; 

END; 

writeln; 

END; 

FOR  m  :=  1  TO  4  DO 
mgradient[m]  :=  -gradient[m]; 

END;  { jacobian  } 

BEGIN  { improve  } 

jacobian(Ecorr,ba,bc,Icorr); 
decomp(order, order J  ac.cond,  ipvt,  work) ; 

IF  (cond  =  cond  +  1) 

THEN 

BEGIN 

IF  detailed  THEN 

writelnC  jacobian  is  singular  to  working  precision  ’); 
FOR  n  :=  1  TO  4  DO 
mgradient[n]  :=  20  *  mgradi  ent[n); 

IF  detailed  THEN 

writelnC  negative  gradient  followed  to  a  new  point '); 
END 
ELSE 
BEGIN 

solve(order,orderJac,mgradient,ipvt); 

END; 

IF  detailed 
THEN 
BEGIN 

writelnC  initial  correction  vector  ’); 

FOR  n  :=  1  TO4DO 
write( mgradient[n]); 
writeln; 

END; 

WHILE  ((abs(mgradient[l])  >  0.5) 

OR  (abs(mgradient[2])  >  0.5) 

OR  (abs(mgradient[3])  >  0.5) 

OR  (abs(mgradient[4])  >  0.5))  DO 


BEGIN 

FOR  i  :=  1  TO  4  DO 
mgradient[i]  :=  mgradient[i]/2; 

IF  detailed  THEN 

write('50%  parameter  change  exceeded1); 

IF  detailed  THEN 

writelnC  reduced  correction  by  half '); 

END;  {  while  } 

newrms  :=  rmserr(Ecorr*(l+mgradient[l]),ba*(l+mgradient[2]),bc*(l+mgradient[3]), 

Icorr*(l+mgradient[4])); 


WHILE  ((newrms  >  rms)  AND  (j  <  25))  DO 
BEGIN 

FOR  i  :=  1  TO  4  DO 
mgradientfi]  :=  mgradientfi]/2; 

IF  detailed  THEN 

writeln('  error  overshoot  reduced  correction  by  half); 
newrms  :=  rmserr(Ecorr*(l+mgradient[l]),ba*(l+mgradient[2]),bc*(l+mgradient[3]), 

Icorr*  ( 1 +mgradient  [4  ] )); 


IF  j  =  25 

THEN  writelnC  no  further  improvement  available '); 

END;  {  while  } 

epsilon  :=  1  +  mgradient[l]; 

alpha  :=  1  +  mgradient[2]; 

beta  :=  1  +  mgradient[3]; 

gamma  :=  1  +  mgradient[4]; 

IF  ((epsilon  =  1)  AND  (alpha  =  1)  AND  (beta  =  1)  AND  (gamma  =  1)) 
THEN 
BEGIN 

writelnC  no  further  improvement  available  '); 

END; 

Ecorr  :=  Ecorr  *  epsilon; 
ba  :=  ba  *  alpha; 
be  :=  be  *  beta; 

Icorr  :=  Icorr  *  gamma; 

IF  ((ba  <=  0)  OR  (be  <=  0)) 

THEN 

BEGIN 

writelnC  negative  tafel  constants  reached '); 


END; 

Rp  :=  ((ba  *  bc)/(ba  +  bc))/(L  *  Icorr); 
writeln; 

writelnC  Ecorr  ba  be  Icorr  Rp'); 
writeln(Ecorr:8;2,ba:  10:2,  be;  10:2,Icorr:  12:6,Rp:12:0); 

END; 

PROCEDURE  stats; 

BEGIN 

j  ;=  0;  m  :=  0;  n  :=  0; 

sdevlcorr  :=  0;  sdevba  :=  0;  sdevbc  :=  0;  sdevRp  :=  0; 

FOR  i  :=  1  TO  count  DO 
BEGIN 

oP  :=  voltdata[i]  -  Ecorr; 

IF  ((abs(oP)  >=  minoP)  AND  (abs(oP)  <=  oPrange)) 

THEN 

BEGIN 

j  >j  +  i; 

IoP  :=  currdata[i); 

tempi  ;=  exp(  L  *  oP/ba); 

temp2  ;=  exp(-L  *  oP/bc); 

sdevlcorr  :=  sdevlcorr  +  sqr(Icorr  -  IoP/(templ  -  temp2)); 

sdevRp  ;=  sdevRp  +  sqr(Rp  -  ((tempi  -  temp2)/(L  *  (1/ba  +  1/bc)  *  IoP))); 


IFoP  >  0 
THEN 
BEGIN 

m  :=  m  +  1; 

sdevba  :=  sdevba  +  sqr(ba  -  L  *  oP/ln(temp2  +  IoP/Icorr)) 
END 
ELSE 

BEGIN 
n  >  n  +  1; 

sdevbc  :=  sdevbc  +  sqr(bc  +  L  *  oP/ln(tempI  -  IoP/Icorr)); 
END; 

END; 

END; 

sdevlcorr  :=  sqrt(sdevIcorr/(j  - 1)); 

sdevRp  :=  sqrt(sdevRp/(j  - 1)); 

sdevba  :=  sqrt(sdevba/(m  - 1)); 

sdevbc  :=  sqrt(sdevbc/(n  -  1)); 

END;  (  stats  } 

BEGIN  { datafit } 


L  :=  ln(10); 

{  This  is  a  Newton's  iteration  loop  for  Ecorr.  The  initial  guess  is  the  rest  potential  Erest.  } 


estEcorr  :=  (potential[l]  +  potential[numofpts])/2; 

REPEAT 

Ecoit  :=  estEcorr; 

estEcorr  :=  Ecorr  -  (eval(Ecorr)/deval(Ecorr)); 

UNTIL  (abs((estEcoir  -  Ecorr)/Ecorr)  <=  IE-6); 

Ecorr  ;=  estEcorr; 

writeln(’  estimated  Ecorr  is  Ecorr  6:2); 

IF  abs(potential[l]  -  Ecorr)  <=  abs(potential[numofpts]  -  Ecorr) 
THEN  oPrange  :=  abs(potential[l]  -  Ecorr) 

ELSE  oPrange  :=  abs(potential[numofpts]  -  Ecor); 
writelnC  range  around  Ecorr  reduced  to  ',oPrange:4:l); 
range  :=oPrange; 
screendata(range); 

suml  ;=  0;  sum2  :=  0;  sum3  :=  0; 
oP  :=  5;  i  :=  0; 

REPEAT 
i  :=  i  +  1; 

Iplus  :=  eval(oP  +  Ecorr); 

Iminus  :=  eval(-oP  +  Ecorr); 
a  :=  ln(abs(IpIus/Iminus))/(2  *  oP); 

IF  odd(i) 

THEN 

BEGIN 

writeC  at  oP  =  ’,oP:5;l); 
write('  a  =  ’,a:8:6); 

END 

ELSE 

BEGIN 

writeC  at  oP  =  ',oP:  5:1); 
writelnC  a  =  \a:8:6); 

END; 

suml  :=  suml  +  a; 
oP  :=  oP  +  0.5; 

UNTIL  oP  >  oPrange; 
num  :=  i; 
a  :=  suml/num; 
writelnC  average  a  is  \a:8:6); 

REPEAT 

write('enter  starting  oP:'); 

writelnC  minimum  oP  where  a  becomes  reasonably  constant'); 
readln;  read(minoP);  ans  :=  99; 

IF  minoP  <=  2 


writelnC  warning  :  do  not  start  at  less  than  2'); 
a ns  :=  98; 

END 

ELSE  IF  minoP  >=  oPrange 
THEN 
BEGIN 

writelnC  warning  :  data  only  good  to  ’,oPrange:4:l); 
ans  :=  98; 

END; 

UNTIL  (ans  <>  98); 

i  :=  0;  suml  ;=  0;  sum2  ;=  0; 
oP  ;=  minoP; 

REPEAT 

Iplus  :=  eval(oP  +  Ecorr); 

Ip3half  :=  eval(1.5*oP  +  Ecorr); 

Iminus  :=  eval(-oP  +  Ecorr); 

Im3half  :=  eval(-1.5*oP  +  Ecorr); 

tempi  :=  sqrt(abs((Ip3half  *  Im3half)/(Iplus  *  Iminus)))/2; 
tempi  :=  (tempi  +  sqrt(sqr(templ)  +  l))/2; 

IF  tempi  >  1.0 
THEN 
BEGIN 

w  :=  2  *  ln(templ  +  sqrt(sqr( tempi)  -  l))/oP; 
i  :=  i  +  1; 

suml  :=  suml  +  w; 

END 

ELSE 

writelnC  data  too  scattered  at  ',oP:4:I,’  to  give  a  w  '); 

oP  :=  oP  +  0.5; 

UNTIL  (oP  >=  2  *  oPrange/3); 
num  ;=  i; 
w  :=  suml/num; 
writelnC  average  w  is  ',w:8:6); 

i  :=  0; 
suml  :=  0; 
oP  :=  minoP; 

REPEAT 

Iplus  :=  eval(oP  +  Ecorr); 

Iminus  :=  eval(-oP  +  Ecorr); 
a  :=  ln(abs(Iplus/Iminus))/(2  *  oP); 

sum2  :=  sum2  +  a; 

Rp  :=  sinh(w  *  oP)/(w  *  (sqri(abs(Iplus  *  Iminus)))); 

IF  detailed  THEN 
suml  :=suml+Rp; 
i  :=  i  +  1; 

oP  :=oP  +  0.5; 

UNTIL  (oP  >=  oPrange); 


a  :=  sum2/num; 

Rp  :=  suml/num; 

writeln('re-estimate  of  a  is  ',a:6:4,'  average  Rp  is  ’*Rp:8:0); 

Icorr  :=  l/(Rp  *  2  *  w); 

writeln('  estimated  Icon  is  '.Icon:  10:8); 

suml  :=  0;  i  :=  0;  oP  :=  -oPrange; 

REPEAT 

IF  (abs(oP)  >=  minoP) 

THEN 

BEGIN 

i  :=  i  +  1; 

IoP  :=  abs(eval(oP  +  Econ)); 

tempi  :=  IoP  1(2  *  Icon  *  exp(a  *  oP)); 
w  :=  (ln(templ  +  sqrt(sqr(  tempi)  +  l)))/(abs(oP)); 

writeln('  at  ',oP:4:l,'  w  =  ',w:8:6); 
suml  :=suml+w; 

END; 
oP  :=  oP  +  0.5; 

UNTIL  (oP  >=  oPrange); 
num  :=  i; 

w  :=  suml/num; 

writeln('re-estimate  of  w  is  ',w:10:8); 

ba  :=  L/(w  +  a);  be  :=  L/(w  -  a);  Icon  :=  1/(2  *  w  *  Rp); 

writeln('  initial  Econ,ba,bc,Icon  and  Rp  are  ’); 

writeln(Econ:6:2,ba:  10:2,bc:  1 0:2, Icon:  12:6, Rp:  1 2:0); 

rms  :=  rmsere(Ecore,ba,bc,Icore); 

writeln('  initial  relative  RMS  enor  is  ',rms:8:6); 

writelnC  enter  1  if  acceptable  '); 

readln;  read(ans); 

WHILE  ans  <>  1  IX) 

BEGIN 

writeln('  enter  ba  be  Rp  *  in  mvolts/decade  and  ohms  *'); 
readln;  read(ba,bc,Rp); 

Icon  :=  ((ba  *  bc)/(ba  +  bc))/(L  *  Rp); 
rms  :=  rmsen(Econ,ba,bc,Icon); 
writeln('relative  RMS  enor  =  ’,rms:8:6); 
writelnC  enter  1  if  acceptable '); 
readln;  read(ans); 

END; 

REPEAT 

writeln('enter  1  to  continue  :  99  to  quit'); 
readln;  read(ans); 

IF  ans  =  1 

THEN  improve(Econ,ba,bc,Icon); 
rms  :=  newrms; 

writeln('relative  RMS  enor  now  =  ',rms:8:6); 


UNTIL  ans  =  99; 
stats; 

writeln('  Ecorr  =  ',Ecorr:6:2,'  millivolts  '); 
writelnO  ba  =  ’,ba:6:2,'  +/-  ',sdevba:6:3,'  millivolts  '); 
writelnC  be  =  ',bc:6:2,'  +/-  ',sdevbc:6:3,'  millivolts  '); 
writeln('  Rp  =  ',Rp:8:0,'  +/-  ’,sdevRp:6: 1,'  ohms  ’); 
writelnC  Icorr  =  ',Icorr:8:6,'  +/-  ',sdevlcornl0:8,'  millamps  '); 

writeln; 

append(LIST); 
writelnC  TAPE  number  ?  '); 
readln;  read(asn); 
writelnC  FILE  number  ? '); 
readln;  read(ans); 
writeln  (LIST); 

writeln(LIST,' - ); 

writeln(LIST,'  File  ’,ans:3); 

writeln(LIST,'  relative  RMS  error  =  ’,rms:8:6); 

writeln(LIST,'  Ecorr  =  ',Ecorr:7:2,'  millivolts  '); 

writeln(LIST,'  Icorr  =  Icorr*  1E3:7:5,'  +/- sdevIcorr*lE3:7:5,'  microamps '); 
writeln(LIST,'  ba  =  ',ba:6:2,'  +/-  ’.sdevbaib^,'  millivolts  ’); 
writeln(LIST,'  be  =  ',bc:6:2,'  +/-  \sdevbc:6:2,'  millivolts  ’); 
writeln(LIST,'  Rp  =  ',Rp*lE-3:6:2,'  +/-  ',sdevRp*lE-3:5:2, '  kohms  ’); 
writeln(LIST,' - ■); 

{  This  section  copies  the  data  and  the  predicted  data  to  the  file  named  PLOTS.,  which  is  acceptable  for  use  } 
{  with  the  plotting  routine  <STAAL>SAPLTF2,  a  simple  curve  plotter.  ) 

rewrite(PLOTS); 
writeln(PLOTS,'CSET  1  O’), 
writeln(PLOTS,'CSET  2  O'); 
writeln(PLOTS,'LSET  1  6'); 
writeln(PLOTS,'LSET  2  O’); 
writeln(PLOTS,'PSET  2.  95.'); 
writeln(PLOTS,'DENSIT  O'); 
writeln(PLOTS, LABEL  1'); 
writeln(PLOTS, 'potential  (millivolts  vs.  SCE)'); 
writeln(PLOTS, LABEL  2'); 

writeln(PLOTS, 'current  density  (microamps/cm!SUP;2!BAK;)'); 
writeln(PLOTS,’LABEL  -4217'); 
writeln(PLOTS,'RMS  error  =  ',rms:6:4); 
writeln(PLOTS;LABEL  -421 1'); 
writeln(PLOTS, 'Ecorr  =  ’,Ecorr:7:2,’  mV  '); 
writeln(PLOTS,'LABEL  -4205'); 

writeln(PLOTS, 'Icorr  =  Icorr*  1E3:7:2,  ’  !FNTl3;m!FNT3;A/cm!SUP;2!BAK; '); 

writeln(PLOTS,'LEGEND’); 

writeln(PLOTS,'data'); 

writeln(PLOTS, 'fitted  curve'); 
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writeln(PLOTS, 'CURVE'); 

FOR  i  :=  1  TO  count  DO 
writeln(PLOTS,voItdata[i],currdata[i]  *  1E3); 
writeln(PLOTS, 'CURVE'); 

FOR  i  :=  1  TO  count  DO 

writeln(PLOTS,voltdata[i],fitcuirent(voltdata[i],Ecorr,ba,bc,Icorr)  *  1E3); 

END; 

BEGIN  {  MAIN  } 
dataheader; 

{  Resets  the  DATA  file  and  reads  in  and  writes  the  information  from  the  start  of  the  file  and  leaves  the  file  } 
{  open  at  the  start  of  the  actual  data.  } 

datagenerator; 

{  Reads  in  the  raw  data  which  is  in  volts  and  namps  and  stores  it  in  arrays  POTENTIAL  and  CURRENT  } 
{  in  mvolts  and  mamps.  POTENTIAL  is  checked  to  ensure  it  is  in  ascending  order,  then  an  interpolating  } 


{  cubic  spline  is  passed  through  the  data.  } 

datafit 

{  Determines  constants  Ecorr,  ba,  be,  Icorr  and  Rp  by  first  making  substitutions  for  the  constants  then  } 

{  separating  the  equation  into  parts  each  with  one  or  two  constants.  Then  solves  each  equation  for  } 

{  approximations  to  its  constant(s).  These  are  then  used  as  the  initial  parameters  to  a  nonlinear  least  } 

{  squares  fitting  routine  which  calculates  the  gradient  and  the  Jacobian  matrix  to  the  sum  of  squared  } 

{  errors  and  generates  a  correction  to  each  parameter.  The  iteration  continues  until  the  size  of  the  } 

{  corrections  becomes  small  enough  that  no  further  improvement  can  be  made.  The  constants  are  then  } 

{  output  with  their  estimated  error  and  the  predicted  data  are  tabulated  with  the  actual  data.  } 


END.  {  main  } 
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4:40pm  Monday,  22  February  1982 


!  <STAAL>C  1 1  A20.CMD 

!  This  sets  up  the  terminal  for  PDP  1 1/34  to  DEC20  connection 

TER  ADM3A 

TER  WIDTH  255 

TER  FORMFEED 

TER  NO  RAISE 

TER  NO  INDICATE 

REFUSE  SYSTEM-MESSAGES 

REFUSE  LINKS 

REFUSE  ADVICE 


C - -DREA-SWA-STAAL-— 

PROGRAM  Cl  1A20 
C 

C  Written  by  Philip  Staal  2:36pm  Thursday,  6  November  1980 
C 

C  PURPOSE:  This  program  types  an  asterisk,  accepts  a  file-name,  types  an  asterisk,  and  accepts  lines  of  ASCII 
C  characters  to  be  put  in  the  file.  The  maximum  line  length  is  255  characters.  When 

C  @@@[CR][LF]  is  received,  the  file  is  closed  and  the  program  stops. 

C 

DOUBLE  PRECISION  FILE 
INTEGER  BUF(255) 

C 

C  ACCEPT  FILE  NAME,  OPEN  FILE,  AND  SEND  * 

C 

TYPE  10010 
ACCEPT  10020, FILE 

OPEN(UNIT=  1  ,DEVICE=DSK',FILE=FILE) 

TYPE  10010 
C 

C  ACCEPT  LINES,  DETERMINE  LENGTH,  AND  WRITE  TO  FILE 
C 

10  ACCEPT  10030, (BUF(J)J=  1,255) 

DO  20  1=1,255 
LAST=256-I 

20  IF(BUF(LAST).NE.32)  GO  TO  30 
30  IF(BUF(LAST).EQ.64)  GO  TO  50 
40  WRITE(  1,1 0030)  (BUF(K),K=1,  LAST) 

GO  TO  10 
C 

C  IF  LAST  CHARACTER  OF  LINE  IS  @,  CHECK  FOR  TWO  PREVIOUS  @'S 


50  IF(BUF(LAST-1).NE.64.0R.BUF(LAST-2).NE.64)  GO  TO  40 
CLOSE(UNIT=  1 ) 

STOP 

C 

10010  FORMATC  *'$) 

10020  FORMAT(A10) 

10030  FORMAT (255R 1 ) 

END 

C - -DREA-SWA-STAAL-- 


•  i 


& 

m 
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This  is  the  format  of  a  typical  data  file  as  it  is  transferred  from  the  EG&G  PARC  Model  350  Corrosion 
Measurement  System  to  the  DREA  DEC  20  computer.  The  name  of  this  file  is  T05F04.DAT.  This  name  also 
refers  to  the  fourth  file  (F)  on  tape  (T)  number  five  with  respect  to  the  tapes  used  in  the  EG&G  PARC  Model  355 
Tape  Storage  Device. 

This  file  is  actually  double  spaced  when  it  is  received  from  the  350  but  is  shown  here  as  single  spaced,  to 
save  space. 

@type  t05f04.dat 

SAMPLE  4 
DATE  28.10 
TAPE  2 
FILE  1 
AREA  1.142E1 
El  -0.276 
EF  -0.227 
MV/SEC  0.167 
ECORR  -0.252 

RESULTS 

RP  5.257E3 
ECORR  -0.256 

-0.276  -5.727E3 
-0.276  -5.063E3 
-0.275  -4.662E3 
-0.275  -4.388E3 
-0.274  -4.046E3 
-0.274  -3.811E3 
-0.273  -3.616E3 
-0.273  -3.440E3 
-0.272  -3.284E3 
-0.272  -3.147E3 
-0.271  -3.020E3 
-0.271  -2.883E3 
-0.270  -2.766E3 
-0.270  -2.658E3 
-0.269 -2.55 1E3 
-0.269  -2.453E3 
-0.268  -2.355E3 
-0.268  -2.267E3 


-0.267  -2.169E3 
-0.267  -2.08 1E3 
-0.266 -1.993E3 
-0.266  -1.905E3 
-0.265  -1.827E3 
-0.265  -1.739E3 
-0.264-1. 65 1E3 
-0.264 -1.573E3 
-0.263  -1.485E3 
-0.263  -1.397E3 
-0.262  -1.299E3 
-0.262  -1.221E3 
-0.261  -1.133E3 
-0.261  -1.036E3 
-0.260  -9.490E2 
-0.260  -8.601E2 
-0.259  -7.633E2 
-0.259  -6.666E2 
-0.258  -5.669E2 
-0.258  -4.662E2 
-0.257  -3.675E2 
-0.257  -2.619E2 
-0.256  -1.612E2 
-0.256  -6.998E1 
-0.255  2.023E1 
-0.255  1.094E2 
-0.254  2.003E2 
-0.254  2.922E2 
-0.253  3.743E2 
-0.253  4.623E2 
-0.252  5.629E2 
-0.252  6.685E2 
-0.251  7.780E2 
-0.251  8.806E2 
-0.250  9.920E2 
-0.250  1.104E3 
-0.249  1.212E3 
-0.249  1.319E3 
-0.248  1.436E3 
-0.248  1.554E3 
-0.247  1.671E3 
-0.247  1.818E3 
-0.246  1.974E3 
-0.246  2.130E3 
-0.245  2.296E3 
-0.245  2.472E3 
-0.244  2.629E3 
-0.244  2.795E3 
-0.243  2.932E3 


-0.243  3.069E3 
-0.242  3.205E3 
-0.242  3.323E3 
-0.241  3.440E3 
-0.241  3.567E3 
-0.240  3.694E3 
-0.240  3.821E3 
-0.239  3.948E3 
-0.239  4.085E3 
-0.238  4.232E3 
-0.238  4.378E3 
-0.237  4.525E3 
-0.237  4.652E3 
-0.236  4.779E3 
-0.236  4.906E3 
-0.235  5.053E3 
-0.235  5.209E3 
-0.234  5.356E3 
-0.234  5.502E3 
-0.233  5.649E3 
-0.233  5.805E3 
-0.232  5.962E3 
-0.232  6.099E3 
-0.231  6.226E3 
-0.231  6.353E3 
-0.230  6.499E3 
-0.230  6.626E3 
-0.229  6.753E3 
-0.229  6.900E3 
-0.228  7.027E3 
-0.228  7.144E3 
-0.227  7.340E3 
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SUBROUTINE  SPLINE  (N,  X,  Y,  B,  C,  D) 
INTEGER  N 

REAL  X(N),  Y(N),  B(N),  C(N),  D(N) 


C 

C  THE  COEFFICIENTS  B(I),  C(I),  AND  D(I),  1=1,2,..., N  ARE  COMPU  TED 
C  FOR  A  CUBIC  INTERPOLATING  SPLINE 
C 

C  S(X)  =  Y(I)  +  B(I)*(X-X(I))  +  C(D*(X-X(I))**2  +  D(I)*(X-X(I))**3 

c 

C  FOR  X(I).LE.X.LE.  X(I+1) 

C 

C  INPUT.. 

C 

C  N  -  THE  NUMBER  OF  DATA  POINTS  OR  KNOTS  (N.GE.2) 

C  X  =  THE  ABSCISSAS  OF  THE  KNOTS  IN  STRICTLY  INCREASING  ORDER 
C  Y  =  THE  ORDINATES  OF  THE  KNOTS 
C 

C  OUTPUT.. 

C 

C  B,  C,  D  =  ARRAYS  OF  SPLINE  COEFFICIENTS  AS  DEFINED  ABOVE. 

C 

C  USING  P  TO  DENOTE  DIFFERENTIATION, 

C 

C  Y(I)  =  S(X(I)) 

C  B(I)  =  SP(X(I)) 

C  C(I)  =  SPP(X(I))/2 

C  D(I)  =  SPPP(X(I))/6  (DERIVATIVE  FROM  THE  RIGHT) 

C 

C  THE  ACCOMPANYING  SUBROUTINE  SEVAL  CAN  BE  USED 
C  TO  EVALUATE  THE  SPLINE  AND  ITS  DERIVATIVE. 

C 

C 

INTEGER  NM1,  IB,  I 
REALT 
C 

NM1  =  N-l 

IF  ( N  .LT.  2  )  RETURN 
IF  ( N  .LT.  3  )  GO  TO  SO 
C 

C  SET  UP  TRIDIAGONAL  SYSTEM 
C 

CB»  DIAGONAL,  D  =  OFFDI  AGONAL,  C  =  RIGHT  HAND  SIDE. 

C 
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D(l)  =  X(2)  -  X(l) 

C(2)  =  (Y(2)  -  Y(l))/D(l) 

DO  10 1  =  2,  NM1 

D(i)  =  x(i+i)-xa) 

B(I)  =  2.*(D(I-1)  +  D(D) 

C(I+1)  =  (Y(I+1)-Y(I))/D(I) 

C(I)  =  C(I+1)  -  C(I) 

10  CONTINUE 
C 

C  END  CONDITIONS.  THIRD  DERIVATIVES  AT  X(l)  AND  X(N) 
C  OBTAINED  FROM  DIVIDED  DIFFERENCES 
C 

B(l)  =  -D(l) 

B(N)  =  -D(N-1) 

C(1)  =  0. 

C(N)  =  0. 

IF  (  N  .EQ.  3  )  GO  TO  15 

C(l)  =  C(3)/(X(4)-X(2))  -  C(2)/(X(3)-X(  1 )) 

C(N)  =  C(N- 1  )/(X(N)-X(N -2))  -  C(N-2)/(X(N-l)-X(N-3)) 

C(l)  =  C(1)*D(1)**2/(X(4)-X(1)) 

C(N)  =  -C(N)*D(N-l)**2/(X(N)-X(N-3)) 

C 

C  FORWARD  ELIMINATION 
C 

15  DO  20 1  =  2,  N 
T  =  D(I-1)/B(I-1) 

B(I)  =  B(I)  -T*D(I-1) 

C(I)-C(I)-T*C(I-1) 

20  CONTINUE 
C 

C  BACK  SUBSTITUTION 
C 

C(N)  =  C(N)/B(N) 

DO  30  IB  =  1,  NM1 
I  =  N-IB 

C(I)  =  (C(I)-D(I)*C(I+1))/B(I) 

30  CONTINUE 
C 

C  COMPUTE  POLYNOMIAL  COEFFICIENTS 
C 

B(N)  =  (Y(N)  -  Y(NM1))/D(NM1)  +  D(NM1)*(C(NM1)  +  2.*C(N)) 
DO  40  I  =  1,  NM1 

B(I)  -  (Y(I+1)  -  Y(I))/D(I)  -  D(I)*(C(I+1)  +  2.*C(I)) 

D(I)  =  (C(I+1)  -  C(I))/D(I) 

C(I)  =  3,*C(I) 

40  CONTINUE 
C(N)  =  3.*C(N) 

D(N)  =  D(N-1) 

RETURN 
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c 

50  B(l)  =  (Y(2)-Y(1))/(X(2)-X(1)) 

C(1)  =  0. 

D(l)  =  0. 

RETURN 

END 

C 

C 

SUBROUTINE  SEVAL(N,  U,  X,  Y,  B,  C,  D.VALDVAL) 

INTEGER  N 

REAL  U,  X(N),  Y(N),  B(N),  C(N),  D(N),  VAL,  DVAL 
C 

C  THIS  SUBROUTINE  EVALUATES  THE  CUBIC  SPLINE  FUNCTION 
C  AND  ITS  DERIVATIVE  AND  RETURNS  VAL  AND  DVAL 
C 

C  VAL  =  Y(I)  +  B(I)*(U-X(I))  +  C(I)*(U-X(I))**2  +  D(I)*(U-X(I))**3 
C  DVAL  =  B(I)  +  2*C(I)*(U-X(I))  +  3*D(1)*(U-X(I))**2 
C  WHERE  X(I)  .LT.  U  .LT.  X(I+1),  USING  HORNER'S  RULE 
C 

C  IF  U  .LT.  X(l)  THEN  I  =  1  IS  USED. 

C  IF  U  .GE.  X(N)  THEN  I  =  N  IS  USED. 

C 

C  INPUT.. 

C 

C  N  =  THE  NUMBER  OF  DATA  POINTS 

C  U  =  THE  ABSCISSA  AT  WHICH  THE  SPLINE  IS  TO  BE  EVALUATED 
C  X,Y  =  THE  ARRAYS  OF  DATA  ABSCISSAS  AND  ORDINATES 
C  B,C,D  =  ARRAYS  OF  SPLINE  COEFFICIENTS  COMPUTED  BY  SPLINE 
C 

C  IF  U  IS  NOT  IN  THE  SAME  INTERVAL  AS  THE  PREVIOUS  CALL,  THEN  A 
C  BINARY  SEARCH  IS  PERFORMED  TO  DETERMINE  THE  PROPER  INTERVAL. 
C 

INTEGER  I,  J,  K 
REAL  DX 
DATA  1/1/ 

IF  ( I  .GE.  N  )  I  =  1 
IF  ( U  .LT.  X(I) )  GO  TO  10 
IF  ( U  .LE.  X(I+1) )  GO  TO  30 
C 

C  BINARY  SEARCH 
C 

10  1  =  1 
J  =  N+l 
20  K  =  (I+J)/2 

IF  (  U  .LT.  X(K) )  J  =  K 
IF  ( U  .GE.  X(K) )  I  =  K 
IF  ( J  .GT.  1+ 1  )  GO  TO  20 
C 

C  EVALUATE  SPLINE  AND  DERIVATIVE 
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30  DX  =  U  -  X(I) 

V  AL  =  Y(I)  +  DX*(B(I)  +  DX*(C(I)  +  DX*D(I))) 

DVAL  =  B(I)  +  DX*(2*C(I)  +  DX*(3*D(I))) 

RETURN 

END 

SUBROUTINE  DECOMP(NDIM,N,A,COND,IPVT,WORK) 

C 

INTEGER  NDIM.N 

REAL  A(NDIM,N),COND,WORK(N) 

INTEGER  IPVT(N) 

C 

C  DECOMPOSES  A  REAL  MATRIX  BY  GAUSSIAN  ELIMINATION 
C  AND  ESTIMATES  THE  CONDITION  OF  THE  MATRIX. 

C 

C  USE  SOLVE  TO  COMPUTE  SOLUTIONS  TO  LINEAR  SYSTEMS. 

C 

C  INPUT.. 

C 

C  NDIM  =  DECLARED  ROW  DIMENSION  OF  THE  ARRAY  CONTAINING  A. 

C 

C  N  =  ORDER  OF  THE  MATRIX. 

C 

C  A  =  MATRIX  TO  BE  TRIANGULARIZED. 

C 

C  OUTPUT.. 

C 

C  A  CONTAINS  AN  UPPER  TRIANGULAR  MATRIX  U  AND  A  PERMUTED 
C  VERSION  OF  A  LOWER  TRIANGULAR  MATRIX  I-L  SO  THAT 
C  (PERMUTATION  MATRIX)* A  =  L*U  . 

C 

C  COND  »  AN  ESTIMATE  OF  THE  CONDITION  OF  A  . 

C  FOR  THE  LINEAR  SYSTEM  A*X  =  B,  CHANGES  IN  A  AND  B 
C  MAY  CAUSE  CHANGES  COND  TIMES  AS  LARGE  IN  X  . 

C  IF  COND  +  1.0  .EQ.  COND,  A  IS  SINGULAR  TO  WORKING  PRECISION. 

C  COND  =  1 .0D+32  IF  EXACT  SINGULARITY  IS  DETECTED. 

C 

C  IPVT  =  THE  PIVOT  VECTOR. 

C  IPVT(K)  =  THE  INDEX  OF  THE  K-TH  PIVOT  ROW 
C  IPVT(N)  =  (-1)**(NUMBER  OF  INTERCHANGES) 

C 

C  WORK  SPACE..  THE  VECTOR  WORK  MUST  BE  DECLARED  AND  INCLUDED 
C  IN  THE  CALL.  ITS  INPUT  CONTENTS  ARE  IGNORED. 

C  ITS  OUTPUT  CONTENTS  ARE  USU  ALL  Y  UN  IMPORTANT. 

C 

C  THE  DETERMINANT  OF  A  CAN  BE  OBTAINED  ON  OUTPUT  BY 
C  DET(A)  =  IPVT(N)  *  A(l,l)  *  A(2,2)  *  ...  *  A(N,N). 

C 

REAL  EK,  T,  ANORM.  YNORM,  ZNORM 
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nnw  non  non  r-  non  o  n  o  -  ,j>  non 


C  REAL  DABS  !ch 

INTEGER  NM1, 1,  J,  K,  KP1,  KB,  KM1,  M 
C 

IPVT(N)  =  1 
EF(NEQ.  1)  GO  TO  80 
NM1  =  N  - 1 

COMPUTE  1-NORM  OF  A 


ANORM  =  0.0D0 
DO  10  J  =  1,N 
T  =  0.0D0 
DO  5  I  =  1,  N 

a(i,j)=dble(a(i,j))  !ch 

T  =  T  +  ABS(A(IJ))  !ch 

CONTINUE 

IF  (T  .GT.  ANORM)  ANORM  =  T 
)  CONTINUE 

GAUSSIAN  ELIMINATION  WITH  PARTIAL  PIVOTING 

DO  35  K  =  1.NM1 
KP1=  K+l 

FIND  PIVOT 

M  =  K 

DO  151  =  KP1,N 

a(i,k)=dble(a(i,k))  !ch 

a(m,k)=dble(a(m,k))  !ch 

IF  (ABS(A(I,K))  .GT.  ABS(A(M,K)))  M  =  I  !ch 
I  CONTINUE 
EPVT(K)  =  M 

IF  (M  .NE.  K)  IPVT(N)  =  -IPVT(N) 

T  =  A(M,K) 

A(M,K)  =  A(K,K) 

A(K,K)  =  T 


SKIP  STEP  IF  PIVOT  IS  ZERO 

IF  (T  .EQ.  0.0D0)  GO  TO  35 

COMPUTE  MULTIPLIERS 

DO  20 1  =  KP1,N 
A(I,K)  =  -A(I,K)/T 
I  CONTINUE 

INTERCHANGE  AND  ELIMINATE  BY  COLUMNS 
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DO  30  J  =  KP1.N 
T  =  A(MJ) 

A(MJ)  =  A(KJ) 

A(K,J)  =  T 

IF  (T  .EQ.  O.ODO)  GO  TO  30 
DO  25  I  =  KP1,N 
A(I,J)  =  A(IJ)  +  A(I,K)*T 
25  CONTINUE 
30  CONTINUE 
35  CONTINUE 
C 

C  COND  =  (1-NORM  OF  A)*(AN  ESTIMATE  OF  1-NORM  OF  A-IN VERSE) 

C  ESTIMATE  OBTAINED  BY  ONE  STEP  OF  INVERSE  ITERATION  FOR  THE 
C  SMALL  SINGULAR  VECTOR.  THIS  INVOLVES  SOLVING  TWO  SYSTEMS 
C  OF  EQUATIONS,  (A-TRANSPOSE)*Y  =  E  AND  A*Z=Y  WHERE  E 
C  IS  A  VECTOR  OF  +1  OR  -1  CHOSEN  TO  CAUSE  GROWTH  IN  Y. 

C  ESTIMATE  =  ( 1  -NORM  OF  Z )/( 1  -NORM  OF  Y) 

C 

C  SOLVE  (A-TRANSPOSE)*Y  =  E 
C 

DO  50  K  =  1,  N 
T  =  O.ODO 

IF  (K  .EQ.  1)  GO  TO  45 
KM1  =  K-l 
DO  40 1  =  L  KM1 
T  =  T  +  A(I,K)*WORK(I) 

40  CONTINUE 
45  EK  =  1 .0D0 

IF  (T  .LT.  O.ODO)  EK  =  -1.0D0 
IF  (A(K,K)  .EQ.  O.ODO)  GO  TO  90 
WORK(K)  =  -(EK  +  T)/A(K,K) 

50  CONTINUE 

DO  60  KB  =  1,  NM1 
K  =  N  -  KB 
T  =  O.ODO 
KP1  =  K+l 
DO  55  I  =  KP1,  N 
T  =  T  +  A(I,K)*WORK(K) 

I  55  CONTINUE 

WORK(K)  =  T 
M  =  IPVT(K) 

IF  (M  .EQ.  K)  GO  TO  60 
T  =  WORK(M) 

WORK(M)  =  WORK(K) 

WORK(K)  =  T 
60  CONTINUE 

YNORM  =  O.ODO 
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DO  65  I  =  1,  N 

work(i)=dble(work(i))  !ch 

YNORM  =  YNORM  +  ABS(WORK(I))  !ch 

65  CONTINUE 

SOLVE  A*Z  =  Y 

CALL  SOLVE (NDIM,  N,  A,  WORK,  EPVT) 

ZNORM  =  0.0D0 
DO  70 1  =  1,  N 

work(i)=dble(work(i))  !ch 

ZNORM  =  ZNORM  +  ABS(WORK(I))  !ch 

)  CONTINUE 

ESTIMATE  CONDITION 

COND  =  ANORM*ZNORM/ YNORM 
IF  (COND  .LT.  1.0D0)  COND  =  1.0D0 
RETURN 


1-BY-l 


)  COND  =  l.ODO 
IF  (A(l,l)  .NE.  0.0D0)  RETURN 

EXACT  SINGULARITY 

)  COND  =  1.0D+32 
RETURN 
END 


SUBROUTINE  SOLVE(NDIM,  N,  A,  B,  IPVT) 

INTEGER  NDIM,  N,  IPVT(N) 

REAL  A(NDIM,N),B(N) 

SOLUTION  OF  LINEAR  SYSTEM,  A*X  =  B  . 

DO  NOT  USE  IF  DECOMP  HAS  DETECTED  SINGULARITY. 

INPUT.. 

NDIM  =  DECLARED  ROW  DIMENSION  OF  ARRAY  CONTAINING  A  . 
N  =  ORDER  OF  MATRIX. 

A  =  TRIANGULARIZED  MATRIX  OBTAINED  FROM  DECOMP . 
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B  =  RIGHT  HAND  SIDE  VECTOR. 

IPVT  =  PIVOT  VECTOR  OBTAINED  FROM  DECOMP 
OUTPUT.. 

B  =  SOLUTION  VECTOR,  X  . 

INTEGER  KB,  KM1,  NM1,  KP1, 1,  K,  M 
REAL  T 

FORWARD  ELIMINATION 

IF  (N  JEQ.  1)  GO  TO  50 
NM1  =  N-l 
EX)  20  K  =  l.NMl 
KP1  =  K+l 
M  =  IPVT(K) 

T  =  B(M) 

B(M)  =  B(K) 

B(K)  =  T 
DO  10 1  =  KP1,  N 
B(I)  =  B(I)  +  A(I,K)*T 
)  CONTINUE 
)  CONTINUE 

BACK  SUBSTITUTION 

DO  40  KB  =  l.NMl 
KM1  =  N-KB 
K  =  KM  1  +  1 
B(K)  =  B(K)/A(K,K) 

T  =  -B(K) 

DO  301=  1,  KM1 
B(I)  =  B(I)  +  A(I,K)*T 
30  CONTINUE 
40  CONTINUE 
50  B(l)  =  B(l)/A(l,l) 

RETURN 


oP,  in  the  following,  is  equivalent  to  A<t>  or  4>  -  <l>corT. 


@ex@corrode 

PASCAL:  CORROS 

FORTRAN:  CORSUB 

SPLINE 

SEVAL 

DECOMP 

SOLVE 

LINK:  Loading 

[LNKXCT  CORROS  execution] 

INPUT  : 

OUTPUT  : 

PLOTS  : 

DATA  :  t05f04.dat 
LIST  : 

Would  you  like  to  read  the  introduction  ?  yes  =  1 
2 

do  you  want  a  detailed  run  ?  yes  =  1 
2 

SAMPLE  4  DATE  28.10  TAPE  2  FILE  1  AREA  1.142E1  El -0. 276  EF -0.227 

MV/SEC  0.167  ECORR -0.252 

RESULTS  RP  5.257E3 

ECORR -0.256  ***  end  of  tape  header  *** 

99  data  points  read  in 

do  you  want  to  see  the  data  tabulated?  yes  =  1 
2 

estimated  Ecorr  is  -255.11 
range  around  Ecorr  reduced  to  20.9 


84  points  in  range 

atoP  =  5.0  a=  -0.000173 

atoP=  5.5 

a  =  0.001862 

at  oP  =  6.0 

a  =  0.002409 

at  oP  =  6.5 

a  =  0.003409 

at  oP  =  7.0 

a  =  0.004673 

atoP=  7.5 

a  =  0.004996 

at  oP  =  8.0 

a  =  0.005475 

atoP=  8.5 

a  =  0.006706 

at  oP  =  9.0 

a  =  0.008307 

atoP=  9.5 

a  =  0.009173 

atoP=  10.0 

a  =  0.010085 

atoP=  10.5 

a  =  0.011198 

at  oP  =  1 1 .0 

a  =  0.011510 

at  oP  =  11.5 

a  =  0.011866 

at  oP  =  12.0 

a  =  0.011721 

at  oP  =  12.5 

a  =  0.011348 

*  *•  *  »  "■"*  •  «  '  » 
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at  oP  =  13.0 

a  =  0.011157 

at  oP  =  13.5 

a  =  0.010626 

atoP=  14.0 

a  =  0.010070 

at  oP  =  14.5 

a  =  0.009549 

atoP=  15.0 

a=  0.009089 

atoP  =  15.5 

a  =  0.008512 

atoP=  16.0 

a  =  0.007840 

atoP=  16.5 

a  =  0.007387 

atoP=  17.0 

a=  0.006931 

atoP  =  17.5 

a  =  0.006360 

atoP  =  18.0 

a=  0.005709 

atoP  =  18.5 

a  =  0.004908 

atoP=  19.0 

a  =  0.003753 

at  oP  =  19.5 

a  =  0.002325 

atoP=  20.0 

a  =  0.001458 

atoP=  20.5 

a  =  -0.000052 

average  a  is  0.00688 1 

enter  starting  oP:  minimum  oP  where  a  becomes  reasonably  constant 
5 

average  w  is  0.057436 

re-estimate  of  a  is  0.0069  average  Rp  is  5217. 
estimated  Icorr  is  0.00166878 
re-estimate  of  w  is  0.05757037 
initial  Ecorr.ba, be, Icorr  and  Rp  are 
-255.11  35.73  45.43  0.001665  5217. 

initial  relative  RMS  error  is  0.082915 
enter  1  if  acceptable 
1 

enter  1  to  continue  :  99  to  quit 
1 

Ecorr  ba  be  Icorr  Rp 

-255.77  42.12  47.78  0.001935  5024. 

relative  RMS  error  now  =  0.073427 
enter  1  to  continue  :  99  to  quit 
1 

Ecorr  ba  be  Icorr  Rp 

-255.86  37.33  37.39  0.001556  5215. 

relative  RMS  error  now  =  0.063821 
enter  1  to  continue  :  99  to  quit 
1 

Ecorr  ba  be  Icorr  Rp 

-256.20  34.92  33.73  0.001376  5415. 

relative  RMS  error  now  =  0.060482 
enter  1  to  continue  :  99  to  quit 
1 

Ecorr  ba  be  Icorr  Rp 

-256.18  34.66  33.74  0.001367 

relative  RMS  error  now  =  0.060409 
enter  1  to  continue  :  99  to  quit 
1 


5430. 


1 


Ecorr  ba  be  Icorr  Rp 

-256.18  34.67  33.74  0.001368  5430. 

relative  RMS  error  now  =  0.060409 
enter  1  to  continue  :  99  to  quit 
1 

Ecorr  ba  be  Icorr  Rp 

-256.18  34.67  33.74  0.001368  5430. 

relative  RMS  error  now  =  0.060409 
enter  1  to  continue  :  99  to  quit 
1 

Ecorr  ba  be  Icorr  Rp 

-256.18  34.67  33.74  0.001368  5430. 

relative  RMS  error  now  =  0.060409 
enter  1  to  continue  :  99  to  quit 
1 

Ecorr  ba  be  Icorr  Rp 

-256.18  34.67  33.74  0.001368  5430. 

relative  RMS  error  now  =  0.060409 
enter  1  to  continue  :  99  to  quit 
1 

no  further  improvement  available 

Ecorr  ba  be  Icorr  Rp 

-256.18  34.67  33.74  0.001368  5430. 

relative  RMS  error  now  =  0.060409 
enter  1  to  continue  :  99  to  quit 
99 

relative  RMS  error  now  =  0.060409 

Ecorr  =  -256.18  millivolts 

ba  =  34.67  +/-  4.640  millivolts 

be  =  33.74  +/-  2.949  millivolts 

Rp  =  5430.  +/-  460.8  ohms 

Icorr  =  0.001368+/-  0.00011417  millamps/cm2 

TAPE  number  ? 

5 

FILE  number  ? 

4 


@start 


INPUT 

OUTPUT 

PLOTS 

DATA 

LIST 


t05f04.dat 


Would  you  like  to  read  the  introduction  ?  yes  =  1 
1 

A  program  to  find  the  corrosion  potential  Ecorr,  the  anodic  and  cathodic  Tafel  constants  ba  and  be,  the  polarization 
resistance  Rp  and  the  corrosion  current  density  Icorr,  from  the  mixed  potential  region  around  Ecorr.  The  data  is  read 
from  polarization  resistance  data  from  the  EG&G  PARC  Model  350  Corrosion  Measurement  System  by  another 
program,  <StaaI>CllA201 2 3 4 5 6 7 8 9 10.  The  program  can  be  run  in  two  modes:  detailed  and  simple.  The  only  difference  being 
the  amount  of  output  to  the  terminal.  In  detailed  mode  the  program  outputs  the  gradient,  the  Jacobian  and  its 
construction,  the  correction  vector  and  its  treatment  at  each  iteration.  The  simple  mode  operates  the  same  routines 
but  suppresses  output 

do  you  want  a  detailed  run  ?  yes  =  1 
1 

SAMPLE  4  DATE  28.10  TAPE  2  FILE  1  AREA  1.142E1  El -0.276  EF -0.227 

MV/SEC  0.167  ECORR -0.252 

RESULTS  RP  5.257E3 

ECORR  -0.256  ***  end  of  tape  header  *** 

99  data  points  read  in 

do  you  want  to  see  the  data  tabulated?  yes  =  1 
1 

************************ ************************ *********** 


Potential  Current  Potential  Current 


1  -276.00  -0.005727 

2  -275.50  -0.005063 

3  -275.00  -0.004662 

4  -274.50  -0.004388 

5  -274.00  -0.004046 

6  -273.50  -0.003811 

7  -273.00  -0.003616 

8  -272.50  -0.003440 

9  -272.00  -0.003284 

10  -271.50  -0.003147 


50  -251.50  0.000669 

51  -251.00  0.000778 

52  -250.50  0.000881 

53  -250.00  0.000992 

54  -249.50  0.001104 

55  -249.00  0.001212 

56  -248.50  0.001319 

57  -248.00  0.001436 

58  -247.50  0.001554 

59  -247.00  0.001671 


11  -271.00  -0.003020  60 

12  -270.50  -0.002883  61 

13  -270.00  -0.002766  62 

14  -269.50  -0.002658  63 

15  -269.00  -0.002551  64 

16  -268.50  -0.002453  65 

17  -268.00  -0.002355  66 

18  -267.50  -0.002267  67 

19  -267.00  -0.002169  68 

20  -266.50  -0.002081  69 

21  -266.00  -0.001993  70 

22  -265.50  -0.001905  71 

23  -265.00  -0.001827  72 

24  -264.50  -0.001739  73 

25  -264.00  -0.001651  74 

26  -263.50  -0.001573  75 

27  -263.00  -0.001485  76 

28  -262.50  -0.001397  77 

29  -262.00  -0.001299  78 

30  -261.50  -0.001221  79 

31  -261.00  -0.001133  80 

32  -260.50  -0.001036  81 

33  -260.00  -0.000949  82 

34  -259.50  -0.000860  83 

35  -259.00  -0.000763  84 

36  -258.50  -0.000667  85 

37  -258.00  -0.000567  86 

38  -257.50  -0.000466  87 

39  -257.0  -0.000368  88 

40  -256.50  -0.000262  89 

41  -256.00  -0.000161  90 

42  -255.50  -0.000070  91 

43  -255.00  0.000020  92 

44  -254.50  0.000109  93 

45  -254.00  0.000200  94 

46  -253.50  0.000292  95 

47  -253.00  0.000374  96 

48  -252.50  0.000462  97 

49  -252.00  0.000563  98 

***************************** 

estimated  Ecorr  is  -255. 1 1 
range  around  Ecorr  reduced  to  20.9 


84  points  in  range 

at  oP  =  5.0 

a  =  -0.000173 

at 

at  oP  =  6.0 

a  -  0.002409 

at 

at  oP  =  7.0 

a  =  0.004673 

at 

at  oP  =  8.0 

a  =  0.005475 

at 

at  oP  =  9.0 

a  =  0.008307 

at 

atoP=  10.0 

a  =  0.010085 

at 

-246.50  0.001818 
-246.00  0.001974 
-245.50  0.002130 
-245.00  0.002296 
-244.50  0.002472 
-244.00  0.002629 
-243.50  0.002795 
-243.00  0.002932 
-242.50  0.003069 
-242.00  0.003205 
-241.50  0.003323 
-241.00  0.003440 
-240.50  0.003567 
-240.00  0.003694 
-239.50  0.003821 
-239.00  0.003948 
-238.50  0.004085 
-238.00  0.004232 
-237.50  0.004378 
-237.00  0.004525 
-236.50  0.004652 
-236.00  0.004779 
-235.50  0.004906 
-235.00  0.005053 
-234.50  0.005209 
-234.00  0.005356 
-233.50  0.005502 
-233.00  0.005649 
-232.50  0.005805 
-232.00  0.005962 
-231.50  0.006099 
-231.00  0.006226 
-230.50  0.006353 
-230.00  0.006499 
-229.50  0.006626 
-229.00  0.006753 
-228.50  0.006900 
-228.00  0.007027 

-227.50  0.007144 

***************************** 


oP=  5.5  a  =  0.001862 

oP=  6.5  a  =  0.003409 

oP  =  7.5  a  =  0.004996 

oP=  8.5  a  =  0.006706 

oP=  9.5  a  =  0.009173 

oP=  10.5  a  =  0.011198 


atoP=  11.0  a  =  0.011510  atoP=  11.5  a=  0.011866 
atoP  =  12.0  a  =  0.011721  atoP=  12.5  a  =  0.011348 

atoP  =  13.0  a  =  0.011157  atoP  =  13.5  a=  0.010626 

atoP  =  14.0  a  =  0.010070  atoP=  14.5  a  =  0.009549 

atoP  =  15.0  a  =  0.009089  atoP=  15.5  a=  0.008512 

atoP=  16.0  a  =  0.007840  atoP  =16.5  a  =  0.007387 

atoP  =  17.0  a  =  0.006931  atoP=  17.5  a=  0.006360 

atoP=  18.0  a  =  0.005709  atoP=  18.5  a=  0.004908 

atoP=  19.0  a  =  0.003753  atoP  =19.5  a  =  0.002325 

at  oP  =  20.0  a  =  0.001458  at  oP  =  20.5  a  =  -0.000052 

average  a  is  0.006881 

enter  starting  oP:  minimum  oP  where  a  becomes  reasonably  constant 
5 

average  w  is  0.057436 

re-estimate  of  a  is  0.0069  average  Rp  is  5217. 
estimated  Icorr  is  0.00166878 
re-estimate  of  w  is  0.05757037 
initial  Ecorr.ba, be, Icorr  and  Rp  are 
-255.11  35.73  45.43  0.001665  5217. 

initial  relative  RMS  error  is  0.082915 
enter  1  if  acceptable 
1 

enter  1  to  continue  :  99  to  quit 
1 

gradient  at  each  parameter  +  del 

5.559737E-04  -2.110495E-05  2.477400E-05  -6.371798E-06 
5.498207E-04  -2.090320E-05  2.471899E-05  -6.459422E-06 
5.500678E-04  -2.091288E-05  2.472174E-05  -6.454742E-06 
5.500761E-04  -2.092529E-05  2.471082E-05  -6.436885E-06 

gradient;  dEcorr,dba,dbc,dJcorr 

5.500043E-04  -2.091514E-05  2.471674E-05  -6.449062E-06 

Jacobian  before  symmetry  check 
5.969370E-01  -1.836815E-02  6.341725E-03  7.173367E-03 
-1.898122E-02  1.193507E-03  2.252136E-04  -1.015496E-03 
5.726088E-03  2.246452E-04  4.996309E-04  -5.922402E-04 
7.726374E-03  -1.036022E-03  -5.680192E-04  1.217643E-03 

Jacobian  matrix 

5.969370E-01  -1.867469E-02  6.033906E-03  7.449870E-03 
-1.867469E-02  1.193507E-03  2.249294E-04  -1.025759E-03 
6.033906E-03  2.249294E-04  4.996309E-04  -5.801297E-04 
7.449870E-03  -1.025759E-03  -5.801297E-04  1.217643E-03 

initial  correction  vector 

5.176613E-03  3.579467E-01  1.037736E-01  3.246053E-01 

error  overshoot  reduced  correction  by  half 
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Ecorr  ba  be  Icorr  Rp 
-255.77  42.12  47.78  0.001935  5024. 

relative  RMS  error  now  =  0.073427 
enter  1  to  continue :  99  to  quit 
1 

gradient  at  each  parameter  +  del 

5.857900E-04  -1.760131E-05  1.349978E-05  3.512548E-06 
5.800501E-04  -1.744107E-05  1.344255E-05  3.444029E-06 
5.802585E-04  -1.744835E-05  1.344500E-05  3.447605E-06 
5.802477E-04  -1.746020E-05  1.343397E-05  3.466087E-06 

gradient:  dEcorr,dba,dbc,dIcorr 

5.801946E-04  -1.745080E-05  1.344010E-05  3.453597E-06 

Jacobian  before  symmetry  check 
5.5953 13E-01  -1.445951E-02  6.390474E-03  5.306356E-03 
-1.505189E-02  9.726819E-04  2.443812E-04  -9.401674E-04 
5.96791  IE-03  2.443358E-04  4.901608E-04  -6.131472E-04 
5.895069E-03  -9.568566E-04  -5.992291E-04  1.248921E-03 

Jacobian  matrix 

5.5953 1 3E  01  -1.475570E-02  6.179192E-03  5.600712E-03 
-1.475570E-02  9.726819E-04  2.443585E-04  -9.485120E-04 
6.179192E-03  2.443585E-04  4.901608E-04  -6.061882E-04 
5.6007 12E-03  -9.485120E-04  -6.061882E-04  1.248921E-03 

initial  correction  vector 

3.298722E-04  -1.137024E-01  -2.174872E-01  -1.961592E-01 

Ecorr  ba  be  Icorr  Rp 
-255.86  37.33  37.39  0.001556  5215. 

relative  RMS  error  now  =  0.063821 
enter  1  to  continue  :  99  to  quit 
1 

gradient  at  each  parameter  +  del 

-4.108566E-04  2.253072E-06  -1.302219E-05  9.335019E-06 
-4.1 7240 IE-04  2.412895E-06  -1 .3 1 3840E-05  9.305I70E-06 
-4.169753E-04  2.404754E-06  -1.313292E-05  9.306781E-06 
-4.170642E-04  2.393229E-06  -1.314816E-05  9.327106E-06 

gradient:  dEcorr,dba,dbc,dIcorr 

-4.170889E-04  2.402509E-06  -1.314064E-05  9.314424E-06 

Jacobian  before  symmetry  check 
6.232323E-01  -1.51 1435E-02  1.136286E-02  2.469824E-03 
-1.494373E-02  1.038634E-03  2.245457E-04  -9.280228E-04 
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1.184525E-02  2.244406E-04  7.718995E-04  -7.513677E-04 
2.059539E-03  -9.254222E-04  -7.642825E-04  1.268177E-03 


Jacobian  matrix 

6.232323E-01  -1.502904E-02  1.160406E-02  2.264682E-03 
-1.502904E-02  1.038634E-03  2.244931E-04  -9.267225E-04 
1.160406E-02  2.24493 IE-04  7.718995E-04  -7.578251E-04 
2.264 68 2E-03  -9.267225E-04  -7.578251E-04  1.263177E-03 

initial  correction  vector 

1.356052E-03  -6.445559E-02  -9.785674E-02  -1.153437E-01 

Ecorr  ba  be  Icorr  Rp 
-256.20  34.92  33.73  0.001376  5415. 

relative  RMS  error  now  =  0.060482 
enter  1  to  continue  :  99  to  quit 
1 

gradient  at  each  parameter  +  del 

-5.223650E-05  1.04327  IE-06  -2.257839E-06  1.009662E-06 
-5.892037E-05  1.223551E-06  -2.384136E-06  9.728229E-07 
-5.862422E-05  1.213895E-06  -2.377776E-06  9.749527E-07 
-5.872406E-05  1.202092E-06  -2.393776E-06  9.949578E-07 

gradient;  dEcorr,dba,dbc,dIcoiT 

-5.875 17  IE-05  1.211821E-06  -2.386212E-06  9.825310E-07 

Jacobian  before  symmetry  check 
6.515212E-01  -1.686581E-02  1.274907E-02  2.765046E-03 
-1.685510E-02  1.172992E-03  2.073932E-04  -9.729263E-04 
1.283736E-02  2.076632E-04  8.436587E-04  -7.563330E-04 
2.713062E-03  -9.708174E-04  -7.578336E-04  1.242671E-03 

Jacobian  matrix 

6.515212E-01  -1.686045E-02  1.279321E-02  2.739054E-03 
-1.686045E-02  1.172992E-03  2.075282E-04  -9.718718E-04 
1 .27932 IE-02  2.075282E-04  8.436587E-04  -7.570833E-04 
2.739054E-03  -9.718718E-04  -7.570833E-04  1.242671E-03 

initial  correction  vector 

-8.115 192E-05  -7.518587E-03  1.824736E-04  -6.380777E-03 

Ecorr  ba  be  Icorr  Rp 
-256.18  34.66  33.74  0.001367  5430. 

relative  RMS  error  now  =  0.060409 
enter  1  to  continue  ;  99  to  quit 
1 


gradient  at  each  parameter  +  del 


7.319809E-06  -2.275095E-07  1.263722E-07  6.491456E-08 
6.366540E-07  -4.393530E-08  2.458801E-09  2.510726E-08 
9.344562E-07  -5.378701E-08  8.758978E-09  2.740973E-08 
8.383090E-07  -6.564670E-08  -7.108129E-09  4.729421E-08 

gradient:  dEcorr,dba,dbc,dIcoiT 

8.084107E-07  -5.585315E-08  3.919034E-10  3.490874E-08 

Jacobian  before  symmetry  check 
6.511399E-01  -1.717567E-02  1.260455E-02  2.989827E-03 
-1.716563E-02  1.191785E-03  2.066145E-04  -9.793553E-04 
1.259803E-02  2.066898E-04  8.367074E-04  -7.500033E-04 
3.000582E-03  -9.801482E-04  -7.49901  IE-04  1.238547E-03 

Jacobian  matrix 

6.511399E-01  -1.717065E-02  1.260129E-02  2.995205E-03 
-1.717065E-02  1.191785E-03  2.066521E-04  -9.797518E-04 
1.260129E-02  2.066521E-04  8.367074E-04  -7.499522E-04 
2.995205E-03  -9.797518E-04  -7.499522E-04  1.238547E-03 

initial  correction  vector 

1.476836E-06  2.615836E-04  1.524016E-04  2.674494E-04 

Ecorr  ba  be  Icorr  Rp 
-256.18  34.67  33.74  0.001368  5430. 

relative  RMS  error  now  =  0.060409 
enter  1  to  continue  :  99  to  quit 
1 

gradient  at  each  parameter  +  del 

6.511249E-06  -1.709626E-07  1.255378E-07  3.002279E-08 

-1.702792E-07  1.248890E-08  1.639860E-09  -9.728382E-09 
1.275062E-07  2.649870E-09  7.948902E-09  -7.441514E-09 
3.123887E-08  -9.21841  IE-09  -7.931849E-09  1.2461 10E-08 

gradient:  dEcorr,dba,dbc,dIcorr 

1.336048E-09  5.812808E-10  -4.255068E-10  6.548362E-11 

Jacobian  before  symmetry  check 
6.5099 13E-01  -1.716153E-02  1.261701E-02  2.990282E-03 
-1.715439E-02  1.190762E-03  2.068589E-04  -9.799692E-04 
1.259633E-02  2.065367E-04  8.374409E-04  -7.506342E-04 
2.99573 IE-03  -9.793865E-04  -7.506998E-04  1.239562E-03 

Jacobian  matrix 

6.5099 13E-01  -1.715796E-02  1.260667E-02  2.993006E-03 
-1.715796E-02  1.190762E-03  2.066978E-04  -9.796779E-04 
1.260667E-02  2.066978E-04  8.374409E-04  -7.506670E-04 
2.993006E-03  -9.796779E-04  -7.506670E-04  1.239562E-03 


initial  correction  vector 

-1.812548E-07  -1.093269E-05  -3.205017E-06  -1.019667E-05 
error  overshoot:  reduced  correction  by  half 
error  overshoot:  reduced  correction  by  half 
error  overshoot:  reduced  correction  by  half 
error  overshoot:  reduced  correction  by  half 
error  overshoot:  reduced  correction  by  half 
error  overshoot:  reduced  correction  by  half 
error  overshoot:  reduced  correction  by  half 

Ecorr  ba  be  Icorr  Rp 
-256.18  34.67  33.74  0.001368  5430. 

relative  RMS  error  now  =  0.060409 
enter  1  to  continue  :  99  to  quit 
1 

gradient  at  each  parameter  +  del 

6.512033E-06  -1.709817E-07  1.255612E-07  3.001824E-08 
-1.694675E-07  1.246798E-08  1.663622E-09  -9.732048E-09 
1.283875E-07  2.625029E-09  7.973966E-09  -7.442935E-09 
3.205787E-08  -9.240864E-09  -7.908795E-09  1.245829E-08 

gradient:  dEcorr,dba,dbc,dIcorr 

2.169145E-09  5.580887E-10  -4.025669E-10  6.369305E-11 

Jacobian  before  symmetry  check 
6.509863E-0I  -1.7I6367E-02  1.262183E-02  2.988872E-03 
-1.715398E-02  1.190989E-03  2.066940E-04  -9.798953E-04 
1.259638E-02  2.066189E-04  8.376533E-04  -7.506229E-04 
2.995455E-03  -9.795741E-04  -7.506628E-04  1.239459E-03 

Jacobian  matrix 

6.509863E-01  -1.715882E-02  1.260910E-02  2.992164E-03 
-1.715882E-02  1.190989E-03  2.066565E-04  -9.797347E-04 
1.260910E-02  2.066565E-04  8.376533E-04  -7.506429E-04 
2.992 164E-03  -9.797347E-04  -7.506429E-04  1.239459E-03 

initial  correction  vector 

-1.811797E-07  -1.083091E-05  -3.162521E-06  -1.009062E-05 
error  overshoot:  reduced  correction  by  half 
error  overshoot  reduced  correction  by  half 

Ecorr  ba  be  Icorr  Rp 
-256.18  34.67  33.74  0.001368  5430. 

relative  RMS  error  now  =  0.060409 
enter  1  to  continue  :  99  to  quit 
1 


gradient  at  each  parameter  +  del 

6.51 1890E-06  -1.711419E-07  1.256683E-07  3.000849E-08 
-1.695817E-07  1.230893E-08  1.772985E-09  -9.744071E-09 
1.282406E-07  2.465868E-09  8.082171E-09  -7.454048E-09 
3.196010E-08  -9.402925E-09  -7.800276E-09  1.244962E-08 

gradient:  dEcorr,dba,dbc,dIcorr 

2.053639E-09  3.975629E-10  -2.939675E-10  5.320544E-11 

Jacobian  before  symmetry  check 
6.509837E-01  1.716353E-02  1.261869E-02  2.990646E-03 

-1.715395E-02  1.191137E-03  2.068305E-04  -9.800488E-04 
1.259622E-02  2.066953E-04  8.376139E-04  -7.506309E-04 
2.995529E-03  -9.797276E-04  -7.507254E-04  1.239641E-03 

Jacobian  matrix 

6.509837E-01  -1.715874E-02  1.260746E-02  2.993087E-03 
-1.715874E-02  1.191137E-03  2.067629E-04  -9.798882E-04 
1.260746E-02  2.067629E-04  8.376139E-04  -7.506781E-04 
2.993087E-03  -9.798882E-04  -7.506781E-04  1.239641E-03 

initial  correction  vector 

-1.353360E-07  -7.883587E-06  -2.179073E-06  -7.267384E-06 

Ecorr  ba  be  Icorr  Rp 
-256.18  34.67  33.74  0.001368  5430. 

relative  RMS  error  now  =  0.060409 
enter  1  to  continue  :  99  to  quit 
1 

gradient  at  each  parameter  +  del 

6.510722E-06  -1.715841E-07  1.259860E-07  2.996003E-08 
-1.708609E-07  1.187419E-08  2.090268E-09  -9.797049E-09 
1.269004E-07  2.031868E-09  8.396494E-09  -7.505719E-09 
3.0648 15E-08  -9.835617E-09  -7.482445E-09  1.239516E-08 

gradient:  dEcorr,dba,dbc,dIcorr 

7.985364E-10  -3.859668E-11  2.371081E-11  1.193712E-12 

Jacobian  before  symmetry  check 
6.509923E-01  -1.716594E-02  1.261019E-02  2.984962E-03 
-1.715455E-02  1.191279E-03  2.070465E-04  -9.797020E-O4 
1.259623E-02  2.066557E-04  8.372783E-04  -7.5061 56E-04 
2.995884E-03  -9.798242E-04  -7.5069 13E-04  1.239397E-03 

Jacobian  matrix 

6.509923E-01  -1.716025E-02  1.26032  IE-02  2.990423E-03 
-1.716025E-02  1.191279E-03  2.06851  IE-04  -9.797631E-04 
1.260321E-02  2.06851  IE-04  8.372783E-04  -7.506534E-04 
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2.990423E-03  -9.79763  IE-04  -7.506534E-04  1.239397E-03 
initial  correction  vector 

2.438097E-09  5.157828E-07  3.653581E-07  6.221716E-07 
enor  overshoot;  reduced  correction  by  half 
error  overshoot;  reduced  correction  by  half 
error  overshoot  reduced  correction  by  half 
error  overshoot  reduced  correction  by  half 
error  overshoot  reduced  correction  by  half 
error  overshoot;  reduced  correction  by  half 
error  overshoot;  reduced  correction  by  half 

no  further  improvement  available 

Ecorr  ba  be  Icorr  Rp 
-256.18  34.67  33.74  0.001368  5430. 

relative  RMS  error  now  =  0.060409 
enter  1  to  continue  ;  99  to  quit 
99 

relative  RMS  error  now  =  0.060409 

Ecorr  =  -256.18  millivolts 

ba  =  34.67  +/-  4.640  millivolts 

be  =  33.74  +/-  2.949  millivolts 

Rp  =  5430.  +/-  460.8  ohms 

Icorr  =  0.001368+/-  0.0001 1417  millamps/cm2 

TAPE  number  ? 

5 

FILE  number  ? 

4 


PLOTS 


This  is  the  format  of  a  PLOTS  file  as  it  is  output  when  running  CORROS.  This  file  will  be  accepted 
when  running  <Staal>SAPLTF^  and  will  result  in  graphic  presentation  of  the  data. 

@type  plots.. 

CSET  1  0 
CSET  2  0 
LSET  1  6 
LSET  2  0 
PSET  2.  95. 

DENSITO 
LABEL  1 

potential  (millivolts  vs.  SCE) 

LABEL  2 

current  density  (microamps/cm!SUP;2!BAK;) 

LABEL -4217 
RMS  error  =  0.0604 
LABEL  -4211 
Ecorr=  -256.18  mV 
LABEL  -4205 

Icorr  =  1.36766  !FNT13;m!FNT3;A/cm!SUP;2!BAK; 

LEGEND 

data 

fitted  curve 
CURVE 

-2.760000E+02  -5.727000E+00 
-2.755000E+02  -5.063000E+00 
-2.750000E+02  -4.662000E+00 
-2.745000E+02  -4.388000E+00 
-2.740000E+02  -4.046000E+00 
-2.735000E+02  -3.811000E+00 
-2.730000E+02  -3.616000E+00 
-2.725000E+02  -3.440000E+00 
-2.720000E+02  -3.284000E+00 
-2.7 1 5000E+02  -3.147000E+00 
-2.710000E+02  -3.020000E+00 
-2.705000E+02  -2.883000E+00 
-2.700000E+02  -2.766000E+00 
-2.695000E+02  -2.658000E+00 
-2.690000E+02  -2.551000E+00 
-2.685000E+02  -2.453000E+00 
-2.680000E+02  -2.355000E+00 
-2.675000E+02  -2.267000E+00 


N) 


-2.670000E+02 
-2.665000E+02 
-2.660000E+02 
-2.655000E+02 
-2.650000E+02 
-2.645000E+02 
-2.640000E+02 
-2.635000E+02 
-2.630000E+02 
-2.625000E+02 
-2.620000E+02 
-2.615000E+02 
-2.610000E+02 
-2.605000E+02 
-2.600000E+02 
•2.595000E+02 
-2.590000E+02 
-2.585000E+02 
-2.580000E+02 
-2.575000E+02 
-2.570000E+02 
-2.565000E+02 
-2.560000E+02 
-2.555000E+02 
-2.550000E+02 
•2.545000E+02 
-2.540000E+02 
-2.535000E+02 
-2.530000E+02 
-2.525000E+02 
.520000E+02 
.515000E+02 
-2.510000E+02 
-2.505000E+02 
-2.500000E+02 
-2.495000E+02 
-2.490000E+02 
-2.485000E+02 
-2.480000E+02 
-2.475000E+02 
-2.470000E+02 
-2.465000E+02 
-2.460000E+02 
-2.455000E+02 
-2.450000E+02 
-2.4450Q0E+02 
-2.440000E+02 
-2.435000E+02 
-2.430000E+02 


-2.169000E+00 

-2.081000E+00 

-1.993000E+00 

-1.905000E+00 

-1.827000E+00 

-1.739000E+00 

-1.651000E+00 

-1.573000E+00 

-1.485000E+00 

-1.397000E+00 

-1.299000E+00 

-1.221000E+00 

-1.133000E+00 

-1.036000E+00 

-9.490000E-01 

-8.601000E-01 

-7.633000E-01 

-6.666000E-01 

-5.669000E-01 

-4.662000E-01 

-3.675000E-01 

-2.619000E-01 

-1.612000E-01 

-6.998000E-02 

2.023000E-02 

1.094000E-01 

2.003000E-01 

2.922000E-01 

3.743000E-01 

4.623000E-01 

5.629000E-01 

6.685000E-01 

7.780000E-01 

8.806000E-01 

9.920000E-01 

1.104000E+00 

1.212000E+00 

1.319000E+00 

1.436000E+00 

1.554000E+00 

1.671000E+00 

1.818000E+00 

1.974000E+00 

2.130000E+00 

2.296000E+00 

2.472000E+00 

2.629000E+00 

2.795000E+00 

2.932000E+00 
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-2.425000E+02 

-2.420000E+02 

-2.415000E+02 

-2.410000E+02 

-2.405000E+02 

-2.400000E+02 

-2.395000E+02 

-2.390000E+02 

-2.385000E+02 

-2.380000E+02 

-2.375000E+02 

-2.370000E+02 

-2.365000E+02 

-2.360000E+02 

-2.355000E+02 

-2.350000E+02 

-2.345000E+02 

CURVE 

-2.760000E+02 

-2.755000E+02 

-2.750000E+02 

-2.745000E+02 

-2.740000E+02 

-2.735000E+02 

-2.730000E+02 

-2.725000E+02 

-2.720000E+02 

-2.715000E+02 

-2.710000E+02 

-2.705000E+02 

-2.700000E+02 

-2.695000E+02 

-2.690000E+02 

-2.685000E+02 

-2.680000E+02 

-2.675000E+02 

-2.670000E+02 

-2.665000E+02 

-2.660000E+02 

-2.655000E+02 

-2.650000E+02 

-2.645000E+02 

-2.640000E+02 

-2.635000E+02 

-2.630000E+02 

-2.625000E+02 

-2.620000E+02 

-2.615000E+02 

-2.610000E+02 


3.069000E+00 

3.205000E+00 

3.323000E+00 

3.440000E+00 

3.567000E+00 

3.694000E+00 

3.821000E+00 

3.948000E+00 

4.085000E+00 

4.232000E+00 

4.378000E+00 

4.525000E+00 

4.652000E+00 

4.779000E+00 

4.906000E+00 

5.053000E+00 

5.209000E+00 

-4.920840E+00 
-4.731086E+00 
-4.546863E+00 
-4.367956E+00 
-4.194163E+00 
-4.025277E+00 
-3.861 105E+00 
-3.701454E+00 
-3.546 143E+00 
-3.394990E+00 
-3.2478 19E+00 
-3.104461E+00 
-2.964752E+00 
-2.828527E+00 
-2.695630E+00 
-2.565907E+00 
-2.439209E+00 
-2.315389E+00 
-2.194 305E+00 
-2.075816E+00 
-1.959787E+00 
-1.846082E+00 
-1.734574E+00 
-1.625.'31E+00 
-1.517630E+00 
-1.41 1944E+00 
-1.307954E+00 
-1. 20554 1E+00 
-1.104587E+00 
-1.004976E+00 
-9.065937E-01 
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-2.605000E+02 

-2.600000E+02 

-2.595000E+02 

-2.590000E+02 

-2.585000E+02 

-2.580000E+02 

-2.575000E+02 

-2.570000E+02 

-2.565000E+02 

-2.560000E+02 

-2.555000E+02 

-2.550000E+02 

-2.545000E+02 

-2.540000E+02 

-2.535000E+02 

-2.530000E+02 

-2.525000E+02 

-2.520000E+02 

-2.515000E+02 

-2.510000E+02 

-2.505000E+02 

-2.500000E+02 

-2.495000E+02 

-2.490000E+02 

-2.485000E+02 

-2.480000E+02 

-2.475000E+02 

-2.470000E+02 

-2.465000E+02 

-2.460000E+02 

-2.455000E+02 

-2.450000E+02 

-2.445000E+02 

-2.440000E+02 

-2.435000E+02 

-2.430000E+02 

-2.425000E+02 

-2.420000E+02 

-2.415000E+02 

-2.410000E+02 

-2.405000E+02 

-2.400000E+02 

-2.395000E+02 

-2.390000E+02 

-2.385000E+02 

-2.380000E+02 

-2.375000E+02 

-2.370000E+02 

-2.365000E+02 


-8.093285E-01 

-7.130686E-01 

-6.177044E-01 

-5.231267E-01 

-4.292280E-01 

-3.35901  IE-01 

-2.430398E-01 

-1.505386E-01 

-5.829243E-02 

3.380389E-02 

1.258544E-01 

2.179627E-01 

3.102341E-01 

4.027723E-01 

4.956825E-01 

5.890696E-01 

6.830383E-01 

7.776944E-01 

8.731452E-01 

9.694982E-01 

1.066861E+00 

1.165343E+00 

1.265055E+00 

1.366108E+00 

1.468617E+00 

1.572696E+00 

I.678460E+00 

1.786029E+00 

1.895523E+00 

2.007064E+00 

2.120775E+00 

2.236786E+00 

2.355223E+00 

2.476220E+00 

2.5999 10E+00 

2.7264 34E+00 

2.855928E+00 

2.988540E+00 

3.124414E+00 

3.263705E+00 

3.406563E+00 

3.553150E+00 

3.703626E+00 

3.858161E+00 

4.016924E+00 

4.180092E+00 

4.347845E+00 

4.520369E+00 

4.697855E+00 
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■2.360000E+02  4.880500E+00 
•2.355000E+02  5.068506E+00 
■2.350000E+02  5.262080E+00 
•2.345000E+02  5.461438E+00 


L1SX 


This  is  the  format  of  a  LIST  file  as  it  is  output  when  running  CORROS.  This  file  is  appended  each  time 
CORROS  is  run  and  may  therefore  serve  as  a  summary  of  results  for  a  group  of  polarization  resistance 
measurements. 

@type  list. 


File  4 

relative  RMS  error  =  0.060409 

Ecorr  =  -256.18  millivolts 

Icorr  =  1.36766  +/-  0.11417  microamps/cm2 

ba  =  34.67  +/-  4.64  millivolts 

be  =  33.74  +/-  2.95  millivolts 

Rp  =  5.43  +/-  0.46  kohms 
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Donahue's35  Data 


@type  tpltst.dat 

this  is  a  file  containing  data  from  a  paper 
on  corrosion  analysis  by  Donahue35 

RESULTS 

ECORR  =  200.0  mV 


0. 

-47.E6 

0.02 

-32.E6 

0.04 

-21. E6 

0.06 

-15.E6 

0.08 

-10.E6 

0.1 

-6.8E6 

0.11 

-5.6E6 

0.12 

-4.6E6 

0.13 

-3.9E6 

0.14 

-3.1E6 

0.15 

-2.6E6 

0.155 

-2.3E6 

0.16 

-2.1E6 

0.165 

-1.8E6 

0.17 

-1.6E6 

0.175 

-1.37E6 

0.18 

-1.15E6 

0.185 

-0.90E6 

0.19 

-0.64E6 

0.195 

-0.35E6 

0.2 

0.00 

0.205 

0.42E6 

0.21 

0.96E6 

0.215 

1.76E6 

0.22 

2.50E6 

0.225 

3.8E6 

0.23 

5.1E6 

0.235 

7.2E6 

0.24 

9.5E6 

0.245 

13.E6 

0.25 

18.E6 

0.26 

32.E6 

0.27 

57.E6 

0.28 

100.E6 

@ex@corrode 

PASCAL:  CORROS 

FORTRAN:  CORSUB 

SPLINE 

SEVAL 

DECOMP 

SOLVE 

LINK:  Loading 

[LNKXCT  CORROS  execution] 
INPUT  : 

OUTPUT  : 

PLOTS  : 

DATA  :  tpltst.dat 
LIST  : 


Would  you  like  to  read  the  introduction  ?  yes  =  1 


do  you  want  a  detailed  run  ?  yes  =  1 


this  is  a  test  file  containing  data  from  a  paper  on  corrosion  analysis  by  Donahue35 
RESULTS 

ECORR  =  200.0  mV  ***  end  of  tape  header  *** 

34  data  points  read  in 

do  you  want  to  see  the  data  tabulated?  yes  =  1 
2 

estimated  Ecorr  is  200.00 
range  around  Ecorr  reduced  to  80.0 
27  points  in  range 


at  oP  =  5.0 

a  = 

0.018232 

atoP-  5.5 

a  = 

0.017987 

at  oP  =  6.0 

a  = 

0.017821 

at  oP  -  6.5 

a  «= 

0.017768 

at  oP  =  7.0 

a  = 

0.017836 

at  oP  =  7.5 

a  = 

0.018021 

at  oP  -  8.0 

a  = 

0.018312 

at  oP  =  8.5 

a  = 

0.018697 

atoP-  9.0 

a  = 

0.019161 

atoP=  9.5 

a  = 

0.019691 

atoP-  10.0 

a  = 

0.020273 

atoP=  10.5 

a  = 

0.020884 

atoP-  11.0 

a  « 

0.021464 

atoP  =  11.5 

a  = 

0.021969 

at  oP  =  12.0 

a  « 

0.022374 

atoP  =  12.5 

a  = 

0.022664 

atoP=  13.0 

a  - 

0.022835 

atoP=  13.5 

a  = 

0.022886 

at  oP  =  14.0 

a  = 

0.022821 

atoP=  14.5 

a  = 

0.022642 

atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
at  oP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  - 
atoP  = 
atoP  = 
atoP  = 
at  oP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
at  oP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  « 
atoP  = 
atoP- 
atoP- 


15.0  a  = 
16.0  a  = 
17.0  a  = 
18.0  a  = 
19.0  a  = 
20.0  a  = 
21.0  a  = 
22.0  a  = 
23.0  a  = 
24.0  a  = 
25.0  a  = 
26.0  a  = 
27.0  a  = 
28.0  a  = 
29.0  a  = 
30.0  a  = 
31.0  a  = 
32.0  a  = 
33.0  a  = 
34.0  a  = 
35.0  a  = 
36.0  a  = 
37.0  a  = 
38.0  a  = 
39.0  a  = 
40.0  a  = 
41.0  a  = 
42.0  a  = 
43.0  a  - 
44.0  a  = 
45.0  a  = 
46.0  a  = 
47.0  a  = 
48.0  a  = 
49.0  a  = 
50.0  a  = 
51.0  a  = 
52.0  a  = 
53.0  a  = 
54.0  a  = 
55.0  a  = 
56.0  a  = 
57.0  a  = 
58.0  a  = 
59.0  a  = 
60.0  a  = 
61.0  a  = 
62.0  a  - 
63.0  a  = 


0.022356 

0.021540 

0.020662 

0.019941 

0.019503 

0.019413 

0.019641 

0.019986 

0.020286 

0.020444 

0.020404 

0.020168 

0.019850 

0.019558 

0.019365 

0.019321 

0.019430 

0.019606 

0.019764 

0.019844 

0.019804 

0.019634 

0.019398 

0.019160 

0.018971 

0.018867 

0.018864 

0.018933 

0.019039 

0.019151 

0.019245 

0.019304 

0.019334 

0.019345 

0.019347 

0.019349 

0.019356 

0.019371 

0.019389 

0.019410 

0.019430 

0.019449 

0.019462 

0.019469 

0.019467 

0.019453 

0.019427 

0.019393 

0.019353 


atoP  =* 
atoP  = 
atoP  = 
atoP  => 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
at  oP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  - 
atoP  = 
at  oP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
at  oP  = 
atoP  « 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP- 
at  oP  = 
atoP  - 
atoP  - 
at  oP  = 
at  oP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
atoP  = 
at  oP  = 
atoP  = 
at  oP  = 
atoP  = 


15.5  a  = 

16.5  a  = 

17.5  a  = 

18.5  a  = 

19.5  a  = 

20.5  a  = 

21.5  a  = 

22.5  a  = 

23.5  a  = 

24.5  a  = 

25.5  a  = 

26.5  a  = 

27.5  a  = 

28.5  a  = 

29.5  a  = 

30.5  a  « 

31.5  a  = 

32.5  a  = 

33.5  a  = 

34.5  a  = 

35.5  a  = 

36.5  a  = 

37.5  a- 

38.5  a  = 

39.5  a  = 

40.5  a  = 

41.5  a  = 

42.5  a- 

43.5  a  = 

44.5  a  = 

45.5  a  = 

46.5  a  = 

47.5  a  = 

48.5  a  = 

49.5  a  = 

50.5  a  = 

51.5  a  = 

52.5  a  = 

53.5  a  = 

54.5  a  = 

55.5  a  = 

56.5  a  = 

57.5  a  = 

58.5  a  = 

59.5  a  = 

60.5  a  = 

61.5  a  = 

62.5  a  = 

63.5  a  = 


0.021975 

0.021093 

0.020272 

0.019681 

0.019413 

0.019499 

0.019811 

0.020149 

0.020387 

0.020451 

0.020304 

0.020012 

0.019696 

0.019446 

0.019322 

0.019361 

0.019516 

0.019692 

0.019817 

0.019841 

0.019732 

0.019520 

0.019276 

0.019057 

0.018907 

0.018854 

0.018892 

0.018983 

0.019096 

0.019202 

0.019279 

0.019322 

0.019341 

0.019347 

0.019347 

0.019352 

0.019363 

0.019379 

0.019399 

0.019420 

0.019440 

0.019456 

0.019467 

0.019469 

0.019461 

0.019441 

0.019411 

0.019373 

0.019333 


atoP=  64.0 

a=  0.019312 

at  oP  =  64.5 

a  =  0.019292 

at  oP  =  65.0 

a=  0.019273 

at  oP  =  65.5 

a=  0.019254 

at  oP  =  66.0 

a  -  0.019236 

at  oP  =  66.5 

a  =  0.019220 

at  oP  =  67.0 

a=  0.019205 

atoP  =  67.5 

a  =  0.019192 

at  oP  =  68.0 

a  =  0.019181 

atoP=  68.5 

a  =  0.019172 

atoP  =  69.0 

a  =  0.019165 

at  oP  =  69.5 

a=  0.019160 

atoP=  70.0 

a  =  0.019158 

at  oP  =  70.5 

a  =  0.019158 

atoP=  71.0 

a  =  0.019160 

atoP  =  71.5 

a  =  0.019164 

at  oP  =  72.0 

a  =  0.019170 

atoP  =  72.5 

a  =  0.019176 

atoP  =  73.0 

a  =  0.019184 

atoP  =  73.5 

a  =  0.019192 

at  oP  =  74.0 

a-  0.019200 

at  oP  =  74.5 

a  =  0.019208 

at  oP  =  75.0 

a  =  0.019216 

atoP  =  75.5 

a  =  0.019224 

atoP=  76.0 

a=  0.019231 

atoP  =  76.5 

a  =  0.019237 

atoP=  77.0 

a  =  0.019243 

at  oP  =  77.5 

a  =  0.019247 

atoP=  78.0 

a  =  0.019250 

atoP  =  78.5 

a  =  0.019251 

atoP=  79.0 

a  =  0.019251 

atoP=  79.5 

a  =  0.019248 

atoP  =  80.0 

a  =  0.019244 

average  a  is  0.019598 

enter  starting  oP:  minimum  oP  where  a  becomes  reasonably  constant 
5 

data  too  scattered  at  5.0  to  give  a  w 
average  w  is  0.039641 

re-estimate  of  a  is  0.0196  average  Rp  is  13. 
estimated  Icorr  is  0.95187600 
re-estimate  of  w  is  0.03980096 
initial  Ecorr,ba,bc,Icorr  and  Rp  are 
200.00  38.76  113.99  0.948041  13. 

initial  relative  RMS  error  is  0.083927 
enter  1  if  acceptable 
1 

enter  1  to  continue  :  99  to  quit 
1 

Ecorr  ba  be  Icorr  Rp 

202.49  37.07  106.57  0.837627  14. 

relative  RMS  error  now  =  0.046343 
enter  1  to  continue  :  99  to  quit 
1 

Ecorr  ba  be  Icorr  Rp 

198.95  40.82  122.23  1.038288  13. 

relative  RMS  error  now  =  0.012931 
enter  1  to  continue :  99  to  quit 
1 

Ecorr  ba  be  Icorr  Rp 

199.38  40.44  121.38  1.025582  13. 

relative  RMS  error  now  =  0.008776 
enter  1  to  continue :  99  to  quit 
1 
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Ecorr  ba  be  Icoit  Rp 
199.47  40.53  122.10  1.034750  13. 

relative  RMS  error  now  =  0.005835 
enter  1  to  continue :  99  to  quit 
1 

Ecorr  ba  be  Icorr  Rp 

200.00  40.57  126.53  1.070331  12. 

relative  RMS  error  now  =  0.0057 16 
enter  1  to  continue :  99  to  quit 
1 

Ecorr  ba  be  Icon  Rp 

199.97  40.57  126.43  1.068972  12. 

relative  RMS  enor  now  =  0.005692 
enter  1  to  continue :  99  to  quit 
1 

Econ  ba  be  Icon  Rp 

200.01  40.57  126.83  1.072053  12. 

relative  RMS  error  now  =  0.005689 
enter  1  to  continue  :  99  to  quit 
1 


Econ  ba  be  Icon  Rp 
200.01  40.57  126.75  1.071465  12. 

relative  RMS  enor  now  =  0.005689 
enter  1  to  continue  :  99  to  quit 
1 

Econ  ba  be  Icon  Rp 

200.00  40.57  126.72  1.071242  12. 

relative  RMS  enor  now  =  0.005689 
enter  1  to  continue  :  99  to  quit 
1 

Econ  ba  be  Icon  Rp 

200.00  40.57  126.70  1.071044  12. 

relative  RMS  enor  now  =  0.005689 
enter  1  to  continue  :  99  to  quit 
1 

no  further  improvement  available 

Econ  ba  be  Icon  Rp 

200.00  40.57  126.70  1.07 1044  12. 

relative  RMS  enor  now  =  0.005689 
enter  1  to  continue  :  99  to  quit 
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relative  RMS  error  now  =  0.005689 

Ecorr  =  200.00  millivolts 

ba  =  40.57+/-  0.712  millivolts 

be  =  126.70+/-  13.147  millivolts 

Rp  =  12.  +/-  0.4  ohms 

Icorr=  1.071044+/-  0.03308610  millamps/cm^ 

TAPE  number  ? 

1 

FILE  number  ? 

1 


APPENDIX  K 


PLOTS  From  CORROS  Analysis  of  Donahue’s35  Data 


@type  plots.. 


CSET  1  0 
CSET  2  0 
LSET  1  6 
LSET2  0 
PSET  2.  95. 

DENSIT  0 
LABEL  1 

potential  (millivolts  vs.  SCE) 

LABEL 2 

current  density  (microamps/cm!SUP;2!BAK;) 

LABEL  -4217 

RMS  error  =  0.0057 

LABEL  4211 

Ecorr  =  200.00  mV 

LABEL  4205 


Icon  =  1071.04400  !FNT13;m!FNT3;A/crn!SUP;2!BAK; 
LEGEND 


data 

fitted  curve 
CURVE 

1.200000E+02  4.600000E+03 
1.300000E+02  -3.900000E+03 
1.400000E+02  -3.100000E+03 
1.500000E+02  -2.600000E+03 
1.550000E+02  -2.300000E+03 
1.600000E+02  -2.100000E+03 


1.650000E+02  -1.800000E+03 
1.700000E+02  -1.600000E+03 
1.750000E+02  -1.370000E+03 
1.800000E+02  -1.150000E+03 
1.850000E+02  -9.000000E+02 


1.900000E+02  -6.400000E+02 


1.950000E+02  -3.500000E+02 
2.000000E+02  0.000000E+00 


2.050000E+02 
2.100000E+02 
2.150000E+02 
2.200Q00E+02 
2.250000E+02 
2.300000E+02 
2.350000EJ  02 


4.200000E+02 

9.600000E+02 

1.760000E+03 

2.500000E+03 

3.800000E+03 

5.100000E+03 

7.200000E+03 
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2.400000E+02 
2.450000E+02 
2.500000E+02 
2.600000E+02 
2.700000E+02 
2.800000E+02 
CURVE 
1.200000E+02 
1.300000E+02 
1.400000E+02 
1 .500000E+02 
1.550000E+02 
1.600000E+02 
1.650000E+02 
1.700000E+02 
1.750000E+02 
1.800000E+02 
1.850000E+02 
1.900000E+02 
1.950000E+02 
2.000000E+02 
2.050000E+02 
2. 100000E+02 
2.150000E+02 
2.200000E+02 
2.250000E+02 
2.300000E+02 
2.350000E+02 
2.400000E+02 
2.450000E+02 
2.500000E+02 
2.600000E+02 
2.700000E+02 
2.800000E+02 


9.500000E+03 

1.300000E+04 

1.800000E+04 

3.200000E+04 

5.700000E+04 

1.000000E+05 

-4.572435E+03 
-3.801960E+03 
-3.151405E+03 
-2.594638E+03 
-2.343237E+03 
-2.105133E+03 
-1.876370E+03 
-1.652412E+03 
-1.427900E+03 
-1.196315E+03 
-9.4955 13E+02 
-6.773440E+02 
-3.665 198E+02 
-9.176887E-04 
4.4451 19E+02 
9.962865E+02 
1.693876E+03 
2.588 196E+03 
3.746600E+03 
5.258308E+03 
7.241592E+03 
9.853335E+03 
1.330172E+04 
1.786308E+04 
3.191243E+04 
5.662904E+04 
1.001740E+05 
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